module dshr_strdata_mod

  ! holds data and methods to advance data models
  ! Obtain the model domain and the stream domain for each stream

  use ESMF             , only : ESMF_Mesh, ESMF_RouteHandle, ESMF_Field, ESMF_FieldBundle
  use ESMF             , only : ESMF_Clock, ESMF_VM, ESMF_VMGet, ESMF_VMGetCurrent
  use ESMF             , only : ESMF_DistGrid, ESMF_SUCCESS, ESMF_MeshGet, ESMF_DistGridGet
  use ESMF             , only : ESMF_VMBroadCast, ESMF_MeshIsCreated, ESMF_MeshCreate
  use ESMF             , only : ESMF_CALKIND_NOLEAP, ESMF_CALKIND_GREGORIAN
  use ESMF             , only : ESMF_CalKind_Flag, ESMF_Time, ESMF_TimeInterval
  use ESMF             , only : ESMF_TimeIntervalGet, ESMF_TYPEKIND_R8, ESMF_FieldCreate
  use ESMF             , only : ESMF_FILEFORMAT_ESMFMESH, ESMF_FieldCreate
  use ESMF             , only : ESMF_FieldBundleCreate, ESMF_MESHLOC_ELEMENT, ESMF_FieldBundleAdd
  use ESMF             , only : ESMF_POLEMETHOD_ALLAVG, ESMF_EXTRAPMETHOD_NEAREST_STOD
  use ESMF             , only : ESMF_REGRIDMETHOD_BILINEAR, ESMF_REGRIDMETHOD_NEAREST_STOD
  use ESMF             , only : ESMF_REGRIDMETHOD_CONSERVE, ESMF_NORMTYPE_FRACAREA, ESMF_NORMTYPE_DSTAREA
  use ESMF             , only : ESMF_ClockGet, operator(-), operator(==), ESMF_CALKIND_NOLEAP
  use ESMF             , only : ESMF_FieldReGridStore, ESMF_FieldRedistStore, ESMF_UNMAPPEDACTION_IGNORE
  use ESMF             , only : ESMF_TERMORDER_SRCSEQ, ESMF_FieldRegrid, ESMF_FieldFill, ESMF_FieldIsCreated
  use ESMF             , only : ESMF_REGION_TOTAL, ESMF_FieldGet, ESMF_TraceRegionExit, ESMF_TraceRegionEnter
  use ESMF             , only : ESMF_LOGMSG_INFO, ESMF_LogWrite
  use shr_kind_mod     , only : r8=>shr_kind_r8, r4=>shr_kind_r4, i2=>shr_kind_I2
  use shr_kind_mod     , only : cs=>shr_kind_cs, cl=>shr_kind_cl, cxx=>shr_kind_cxx
  use shr_sys_mod      , only : shr_sys_abort
  use shr_const_mod    , only : shr_const_pi, shr_const_cDay, shr_const_spval
  use shr_cal_mod      , only : shr_cal_calendarname, shr_cal_timeSet
  use shr_cal_mod      , only : shr_cal_noleap, shr_cal_gregorian
  use shr_cal_mod      , only : shr_cal_date2ymd, shr_cal_ymd2date, shr_cal_leapyear
  use shr_orb_mod      , only : shr_orb_decl, shr_orb_cosz, shr_orb_undef_real
#ifdef CESMCOUPLED
  use shr_pio_mod      , only : shr_pio_getiosys, shr_pio_getiotype, shr_pio_getioformat
#endif
  use shr_string_mod   , only : shr_string_listgetname, shr_string_listisvalid, shr_string_listgetnum

  use dshr_stream_mod  , only : shr_stream_streamtype, shr_stream_getModelFieldList, shr_stream_getStreamFieldList
  use dshr_stream_mod  , only : shr_stream_taxis_cycle, shr_stream_taxis_extend, shr_stream_findBounds
  use dshr_stream_mod  , only : shr_stream_getCurrFile, shr_stream_setCurrFile, shr_stream_getMeshFilename
  use dshr_stream_mod  , only : shr_stream_init_from_inline, shr_stream_init_from_esmfconfig
#ifndef DISABLE_FoX
  use dshr_stream_mod  , only : shr_stream_init_from_xml
#endif
  use dshr_stream_mod  , only : shr_stream_getnextfilename, shr_stream_getprevfilename, shr_stream_getData
  use dshr_tinterp_mod , only : shr_tInterp_getCosz, shr_tInterp_getAvgCosz, shr_tInterp_getFactors
  use dshr_methods_mod , only : dshr_fldbun_getfldptr, dshr_fldbun_getfieldN, dshr_fldbun_fldchk, chkerr
  use dshr_methods_mod , only : dshr_fldbun_diagnose, dshr_fldbun_regrid, dshr_field_getfldptr

  use pio              , only : file_desc_t, iosystem_desc_t, io_desc_t, var_desc_t
  use pio              , only : pio_openfile, pio_closefile, pio_nowrite
  use pio              , only : pio_seterrorhandling, pio_initdecomp, pio_freedecomp
  use pio              , only : pio_inquire, pio_inq_varid, pio_inq_varndims, pio_inq_vardimid
  use pio              , only : pio_inq_dimlen, pio_inq_vartype, pio_inq_dimname, pio_inq_dimid
  use pio              , only : pio_double, pio_real, pio_int, pio_offset_kind, pio_get_var
  use pio              , only : pio_read_darray, pio_setframe, pio_fill_double, pio_get_att
  use pio              , only : PIO_BCAST_ERROR, PIO_RETURN_ERROR, PIO_NOERR, PIO_INTERNAL_ERROR, PIO_SHORT

  implicit none
  private

  public  :: shr_strdata_type
  public  :: shr_strdata_init_from_config
  public  :: shr_strdata_init_from_inline
  public  :: shr_strdata_setOrbs
  public  :: shr_strdata_advance
  public  :: shr_strdata_get_stream_domain  ! public since needed by dshr_mod
  public  :: shr_strdata_get_stream_pointer ! get a pointer into a stream's fldbun_model field bundle
  public  :: shr_strdata_get_stream_count
  public  :: shr_strdata_get_stream_fieldbundle
  public  :: shr_strdata_print

  private :: shr_strdata_init_model_domain
  private :: shr_strdata_get_stream_nlev
  private :: shr_strdata_readLBUB

  interface shr_strdata_get_stream_pointer
     module procedure shr_strdata_get_stream_pointer_1d
     module procedure shr_strdata_get_stream_pointer_2d
  end interface shr_strdata_get_stream_pointer

  ! public data members:
  integer                              :: debug = 0  ! local debug flag
  character(len=*) ,parameter, public  :: shr_strdata_nullstr = 'null'
  character(len=*) ,parameter          :: shr_strdata_unset = 'NOT_SET'
  integer          ,parameter          :: main_task = 0

  ! note that the fields in fldbun_stream_lb and fldbun_stream_ub contain the the names fldlist_model

  type shr_strdata_perstream
     character(CL)                       :: stream_meshfile                 ! stream mesh file from stream txt file
     type(ESMF_Mesh)                     :: stream_mesh                     ! stream mesh created from stream mesh file
     type(io_desc_t)                     :: stream_pio_iodesc               ! stream pio descriptor
     logical                             :: stream_pio_iodesc_set =.false.  ! true=>pio iodesc has been set
     type(ESMF_RouteHandle)              :: routehandle                     ! stream n -> model mesh mapping
     character(CL), allocatable          :: fldlist_stream(:)               ! names of stream file fields
     character(CL), allocatable          :: fldlist_model(:)                ! names of stream model fields
     integer                             :: stream_nlev                     ! number of vertical levels in stream
     integer                             :: stream_lb                       ! index of the Lowerbound (LB) in fldlist_stream
     integer                             :: stream_ub                       ! index of the Upperbound (UB) in fldlist_stream
     type(ESMF_Field)                    :: field_stream                    ! a field on the stream data domain
     type(ESMF_Field)                    :: field_stream_vector             ! a vector field on the stream data domain
     type(ESMF_FieldBundle), allocatable :: fldbun_data(:)                  ! stream field bundle interpolated to model grid spatially
     type(ESMF_FieldBundle)              :: fldbun_model                    ! stream n field bundle interpolated to model grid and time
     integer                             :: ymdLB = -1                      ! stream ymd lower bound
     integer                             :: todLB = -1                      ! stream tod lower bound
     integer                             :: ymdUB = -1                      ! stream ymd upper bound
     integer                             :: todUB = -1                      ! stream tod upper bound
     real(r8)                            :: dtmin = 1.0e30_r8
     real(r8)                            :: dtmax = 0.0_r8
     logical                             :: override_annual_cycle = .false.
     type(ESMF_Field)                    :: field_coszen                    ! needed for coszen time interp
  end type shr_strdata_perstream

  type shr_strdata_type
     type(shr_strdata_perstream), allocatable :: pstrm(:)              ! stream info
     type(shr_stream_streamType), pointer :: stream(:)=> null()        ! stream datatype
     logical                        :: mainproc
     integer                        :: io_type                         ! pio info
     integer                        :: io_format                       ! pio info
     integer                        :: modeldt = 0                     ! model dt in seconds
     type(ESMF_Mesh)                :: model_mesh                      ! model mesh
     real(r8), pointer              :: model_lon(:) => null()          ! model longitudes
     real(r8), pointer              :: model_lat(:) => null()          ! model latitudes
     integer                        :: model_nxg                       ! model global domain lon size
     integer                        :: model_nyg                       ! model global domain lat size
     integer                        :: model_nzg                       ! model global domain vertical size
     integer                        :: model_lsize                     ! model local domain size
     integer, pointer               :: model_gindex(:)                 ! model global index spzce
     integer                        :: model_gsize                     ! model global domain size
     type(ESMF_CLock)               :: model_clock                     ! model clock
     character(CL)                  :: model_calendar = shr_cal_noleap ! model calendar for ymd,tod
     integer                        :: ymd, tod                        ! model time
     type(iosystem_desc_t), pointer :: pio_subsystem => null()         ! pio info
     real(r8)                       :: eccen  = SHR_ORB_UNDEF_REAL     ! cosz t-interp info
     real(r8)                       :: mvelpp = SHR_ORB_UNDEF_REAL     ! cosz t-interp info
     real(r8)                       :: lambm0 = SHR_ORB_UNDEF_REAL     ! cosz t-interp info
     real(r8)                       :: obliqr = SHR_ORB_UNDEF_REAL     ! cosz t-interp info
     real(r8), allocatable          :: tavCoszen(:)                    ! cosz t-interp data
  end type shr_strdata_type

  real(r8)         ,parameter :: deg2rad = SHR_CONST_PI/180.0_r8
  character(*)     ,parameter :: u_FILE_u = &
       __FILE__

!===============================================================================
contains
!===============================================================================

  integer function shr_strdata_get_stream_count(sdat)
    type(shr_strdata_type)     , intent(in) :: sdat
    if(associated(sdat%stream)) then
       shr_strdata_get_stream_count = size(sdat%stream)
    else
       shr_strdata_get_stream_count = 0
    endif
  end function shr_strdata_get_stream_count

  !===============================================================================
  type(ESMF_FieldBundle) function shr_strdata_get_stream_fieldbundle(sdat, ns, name)

    ! input/output variables
    type(shr_strdata_type)     , intent(in) :: sdat
    integer                    , intent(in) :: ns ! stream number
    character(len=*)           , intent(in) :: name

    if(trim(name) .eq. 'model') then
       shr_strdata_get_stream_fieldbundle = sdat%pstrm(ns)%fldbun_model
    else if (trim(name) .eq. 'model_lb') then
       shr_strdata_get_stream_fieldbundle = sdat%pstrm(ns)%fldbun_data(sdat%pstrm(ns)%stream_lb)
    else if (trim(name) .eq. 'model_ub') then
       shr_strdata_get_stream_fieldbundle = sdat%pstrm(ns)%fldbun_data(sdat%pstrm(ns)%stream_ub)
    else if (trim(name) .eq. 'stream_lb') then
       call shr_sys_abort("should not be here")
    else if (trim(name) .eq. 'stream_ub') then
       call shr_sys_abort("should not be here")
    else
       call shr_sys_abort(trim(name)//' is not a recognized stream bundle name')
    endif

  end function shr_strdata_get_stream_fieldbundle

  !===============================================================================
  subroutine shr_strdata_init_from_config(sdat, streamfilename, model_mesh, clock, compname, logunit, rc)

    ! input/output variables
    type(shr_strdata_type)     , intent(inout) :: sdat
    character(len=*)           , intent(in)    :: streamfilename
    type(ESMF_Mesh)            , intent(in)    :: model_mesh
    type(ESMF_Clock)           , intent(in)    :: clock
    character(len=*)           , intent(in)    :: compname
    integer                    , intent(in)    :: logunit
    integer                    , intent(out)   :: rc

    ! local variables
    type(ESMF_VM) :: vm
    integer       :: localPet
    character(len=*), parameter  :: subname='(shr_strdata_init_from_config)'
    ! ----------------------------------------------
    rc = ESMF_SUCCESS
    call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO)

#ifdef CESMCOUPLED
    ! Initialize sdat  pio
    sdat%pio_subsystem => shr_pio_getiosys(trim(compname))
    sdat%io_type       =  shr_pio_getiotype(trim(compname))
    sdat%io_format     =  shr_pio_getioformat(trim(compname))
#endif

    call ESMF_VMGetCurrent(vm, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    call ESMF_VMGet(vm, localPet=localPet, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

    ! Initialize sdat streams (read xml file for streams)
    sdat%mainproc = (localPet == main_task)

#ifdef DISABLE_FoX
    call shr_stream_init_from_esmfconfig(streamfilename, sdat%stream, logunit, &
         sdat%pio_subsystem, sdat%io_type, sdat%io_format, rc=rc)
#else
    call shr_stream_init_from_xml(streamfilename, sdat%stream, sdat%mainproc, logunit, &
         sdat%pio_subsystem, sdat%io_type, sdat%io_format, trim(compname), rc=rc)
#endif

    allocate(sdat%pstrm(shr_strdata_get_stream_count(sdat)))

    ! Initialize sdat model domain
    sdat%model_mesh = model_mesh
    call shr_strdata_init_model_domain(sdat, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

    ! Now finish initializing sdat
    call shr_strdata_init(sdat, clock, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

  end subroutine shr_strdata_init_from_config

  !===============================================================================
  subroutine shr_strdata_init_from_inline(sdat, my_task, logunit, compname, &
       model_clock, model_mesh, &
       stream_meshfile, stream_lev_dimname, stream_mapalgo, &
       stream_filenames, stream_fldlistFile, stream_fldListModel, &
       stream_yearFirst, stream_yearLast, stream_yearAlign, &
       stream_offset, stream_taxmode, stream_dtlimit, stream_tintalgo, stream_name, rc)

    ! input/output variables
    type(shr_strdata_type) , intent(inout) :: sdat                   ! stream data type
    integer                , intent(in)    :: my_task                ! my mpi task
    integer                , intent(in)    :: logunit                ! stdout logunit
    character(len=*)       , intent(in)    :: compname               ! component name (e.g. ATM, OCN, ...)
    type(ESMF_Clock)       , intent(in)    :: model_clock            ! model clock
    type(ESMF_Mesh)        , intent(in)    :: model_mesh             ! model mesh
    character(*)           , intent(in)    :: stream_meshFile        ! full pathname to stream mesh file
    character(*)           , intent(in)    :: stream_lev_dimname     ! name of vertical dimension in stream
    character(*)           , intent(in)    :: stream_mapalgo         ! stream mesh -> model mesh mapping type
    character(*)           , intent(in)    :: stream_filenames(:)    ! stream data filenames (full pathnamesa)
    character(*)           , intent(in)    :: stream_fldListFile(:)  ! file field names, colon delim list
    character(*)           , intent(in)    :: stream_fldListModel(:) ! model field names, colon delim list
    integer                , intent(in)    :: stream_yearFirst       ! first year to use
    integer                , intent(in)    :: stream_yearLast        ! last  year to use
    integer                , intent(in)    :: stream_yearAlign       ! align yearFirst with this model year
    integer                , intent(in)    :: stream_offset          ! offset in seconds of stream data
    character(*)           , intent(in)    :: stream_taxMode         ! time axis mode
    real(r8)               , intent(in)    :: stream_dtlimit         ! ratio of max/min stream delta times
    character(*)           , intent(in)    :: stream_tintalgo        ! time interpolation algorithm
    character(*), optional , intent(in)    :: stream_name            ! name of stream
    integer                , intent(out)   :: rc                     ! error code
    ! ----------------------------------------------

    rc = ESMF_SUCCESS

    ! Initialize sdat%logunit and sdat%mainproc
    sdat%mainproc = (my_task == main_task)
#ifdef CESMCOUPLED
    ! Initialize sdat pio
    sdat%pio_subsystem => shr_pio_getiosys(trim(compname))
    sdat%io_type       =  shr_pio_getiotype(trim(compname))
    sdat%io_format     =  shr_pio_getioformat(trim(compname))
#endif

    ! Initialize sdat%pstrm - ASSUME only 1 stream
    allocate(sdat%pstrm(1))

    ! Initialize sdat model domain
    sdat%model_mesh = model_mesh
    call shr_strdata_init_model_domain(sdat, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

    ! Initialize sdat stream - ASSUME only 1 stream
    call shr_stream_init_from_inline(sdat%stream, &
         sdat%pio_subsystem, sdat%io_type, sdat%io_format, &
         stream_meshfile, stream_lev_dimname, stream_mapalgo, &
         stream_yearFirst, stream_yearLast, stream_yearAlign, &
         stream_offset, stream_taxmode, stream_tintalgo, stream_dtlimit, &
         stream_fldlistFile, stream_fldListModel, stream_fileNames, &
         logunit, trim(compname))

    ! Now finish initializing sdat
    call shr_strdata_init(sdat, model_clock, stream_name, rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

  end subroutine shr_strdata_init_from_inline

  !===============================================================================
  subroutine shr_strdata_init_model_domain( sdat, rc)

    ! ----------------------------------------------
    ! Initialize sdat model domain info
    ! ----------------------------------------------

    ! input/output variables
    type(shr_strdata_type)     , intent(inout) :: sdat
    integer                    , intent(out)   :: rc

    ! local variables
    integer               :: n                  ! generic counters
    type(ESMF_DistGrid)   :: distGrid
    integer               :: tileCount
    integer, allocatable  :: elementCountPTile(:)
    integer               :: spatialDim         ! number of dimension in mesh
    integer               :: numOwnedElements   ! local size of mesh
    real(r8), allocatable :: ownedElemCoords(:) ! mesh lat and lons
    character(len=*), parameter  :: subname='(shr_strdata_init_model_domain)'
    ! ----------------------------------------------

    rc = ESMF_SUCCESS

    ! initialize sdat%lsize
    call ESMF_MeshGet(sdat%model_mesh, elementdistGrid=distGrid, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    call ESMF_DistGridGet(distGrid, localDe=0, elementCount=sdat%model_lsize, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

    ! initialize sdat%model_gindex
    allocate(sdat%model_gindex(sdat%model_lsize))
    call ESMF_DistGridGet(distGrid, localDe=0, seqIndexList=sdat%model_gindex, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

    ! initialize sdat%model_gsize
    call ESMF_DistGridGet(distGrid, tileCount=tileCount, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    allocate(elementCountPTile(tileCount))
    call ESMF_distGridGet(distGrid, elementCountPTile=elementCountPTile, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    sdat%model_gsize = 0
    do n = 1,size(elementCountPTile)
       sdat%model_gsize = sdat%model_gsize + elementCountPTile(n)
    end do
    deallocate(elementCountPTile)

    ! determine sdat%model_lon, sdat%model_lat
    call ESMF_MeshGet(sdat%model_mesh, spatialDim=spatialDim, &
         numOwnedElements=numOwnedElements, elementdistGrid=distGrid, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    allocate(ownedElemCoords(spatialDim*numOwnedElements))
    call ESMF_MeshGet(sdat%model_mesh, ownedElemCoords=ownedElemCoords)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    allocate(sdat%model_lon(numOwnedElements))
    allocate(sdat%model_lat(numOwnedElements))
    do n = 1, numOwnedElements
       sdat%model_lon(n) = ownedElemCoords(2*n-1)
       sdat%model_lat(n) = ownedElemCoords(2*n)
    end do

  end subroutine shr_strdata_init_model_domain

  !===============================================================================
  subroutine shr_strdata_init(sdat, model_clock, stream_name, rc)

    ! input/output variables
    type(shr_strdata_type) , intent(inout), target :: sdat
    type(ESMF_Clock)       , intent(in)    :: model_clock
    character(*), optional , intent(in)    :: stream_name
    integer                , intent(out)   :: rc

    ! local variables
    type(ESMF_Mesh), pointer     :: stream_mesh
    type(ESMF_CalKind_Flag)      :: esmf_caltype    ! esmf calendar type
    character(CS)                :: calendar        ! calendar name
    integer                      :: ns              ! stream index
    integer                      :: m               ! generic index
    character(CL)                :: fileName        ! generic file name
    integer                      :: nfld            ! loop stream field index
    type(ESMF_Field)             :: lfield          ! temporary
    type(ESMF_Field)             :: lfield_dst      ! temporary
    integer                      :: srcTermProcessing_Value = 0 ! should this be a module variable?
    integer                      :: localpet
    logical                      :: fileExists
    type(ESMF_VM)                :: vm
    logical                      :: mainproc
    integer                      :: nvars
    integer                      :: i, stream_nlev, index
    character(CL)                :: stream_vector_names
    character(len=*), parameter  :: subname='(shr_sdat_init)'
    ! ----------------------------------------------

    rc = ESMF_SUCCESS

    call ESMF_VmGetCurrent(vm, rc=rc)
    call ESMF_VMGet(vm, localpet=localPet, rc=rc)
    mainproc= (localPet==main_task)

    ! Loop over streams
    do ns = 1,shr_strdata_get_stream_count(sdat)

       ! Initialize calendar for stream n
       call ESMF_VMBroadCast(vm, sdat%stream(ns)%calendar, CS, 0, rc=rc)

       ! Set pointer for stream_mesh
       stream_mesh => sdat%pstrm(ns)%stream_mesh

       ! Create the target stream mesh from the stream mesh file
       ! TODO: add functionality if the stream mesh needs to be created from a grid
       call shr_stream_getMeshFileName (sdat%stream(ns), filename)
       if (filename /= 'none' .and. mainproc) then
          inquire(file=trim(filename),exist=fileExists)
          if (.not. fileExists) then
             write(sdat%stream(1)%logunit,'(a)') "ERROR: file does not exist: "//trim(fileName)
             call shr_sys_abort(subName//"ERROR: file does not exist: "//trim(fileName))
          end if
       endif
       !
       ! We do not yet have mask information, but we are required to set it here and change it
       ! later.
       !
       if(filename /= 'none') then
          stream_mesh = ESMF_MeshCreate(trim(filename), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc)
          if (ChkErr(rc,__LINE__,u_FILE_u)) return
       endif

       ! Determine the number of stream levels
       call shr_strdata_get_stream_nlev(sdat, ns, rc=rc)
       stream_nlev = sdat%pstrm(ns)%stream_nlev

       ! Determine field names for stream fields with both stream file names and data model names
       nvars = sdat%stream(ns)%nvars

       ! Allocate memory
       allocate(sdat%pstrm(ns)%fldList_model(nvars))
       call shr_stream_getModelFieldList(sdat%stream(ns), sdat%pstrm(ns)%fldlist_model)
       allocate(sdat%pstrm(ns)%fldlist_stream(nvars))
       call shr_stream_getStreamFieldList(sdat%stream(ns), sdat%pstrm(ns)%fldlist_stream)

       ! Create field bundles on model mesh
       if (sdat%stream(ns)%readmode=='single') then
          sdat%pstrm(ns)%stream_lb = 1
          sdat%pstrm(ns)%stream_ub = 2
          allocate(sdat%pstrm(ns)%fldbun_data(2))
          if (mainproc) then
             write(sdat%stream(1)%logunit,'(a,i8)') trim(subname)//" Creating field bundle array fldbun_data of size 2 for stream ",&
                  ns
          end if
       else if(sdat%stream(ns)%readmode=='full_file') then
          ! TODO: add this in
       endif

       ! Create spatially interpolated (but not time interpolated) field bundle - fldbun_data array
       do i=1,size(sdat%pstrm(ns)%fldbun_data)
          sdat%pstrm(ns)%fldbun_data(i) = ESMF_FieldBundleCreate(rc=rc) ! stream mesh
       enddo
       do nfld = 1, nvars
          do i=1,size(sdat%pstrm(ns)%fldbun_data)
             if (sdat%pstrm(ns)%stream_nlev > 1) then
                lfield = ESMF_FieldCreate(sdat%model_mesh, ESMF_TYPEKIND_r8, &
                     name=trim(sdat%pstrm(ns)%fldlist_model(nfld)), &
                     ungriddedLbound=(/1/), ungriddedUbound=(/stream_nlev/), gridToFieldMap=(/2/), &
                     meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
             else
                lfield = ESMF_FieldCreate(sdat%model_mesh, ESMF_TYPEKIND_r8, &
                     name=trim(sdat%pstrm(ns)%fldlist_model(nfld)), &
                     meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
             end if
             call ESMF_FieldBundleAdd(sdat%pstrm(ns)%fldbun_data(i), (/lfield/), rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
             if (mainproc) then
                if (i == 1) then
                   write(sdat%stream(1)%logunit,'(a,i8)') "    adding field "//trim(sdat%pstrm(ns)%fldlist_model(nfld))//&
                        " to fldbun_data for stream ",ns
                end if
             end if
          enddo
       end do

       ! Create spatially and time interpolated field bundle - fldbun_model
       sdat%pstrm(ns)%fldbun_model = ESMF_FieldBundleCreate(rc=rc) ! time interpolation on model mesh
       do nfld = 1, nvars
          ! create temporary fields on model mesh and add the fields to the field bundle
          if (stream_nlev > 1) then
             lfield = ESMF_FieldCreate(sdat%model_mesh, ESMF_TYPEKIND_r8, &
                  ungriddedLbound=(/1/), ungriddedUbound=(/stream_nlev/), gridToFieldMap=(/2/), &
                  name=trim(sdat%pstrm(ns)%fldlist_model(nfld)), &
                  meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
          else
             lfield = ESMF_FieldCreate(sdat%model_mesh, ESMF_TYPEKIND_r8, &
                  name=trim(sdat%pstrm(ns)%fldlist_model(nfld)), &
                  meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
          end if
          call ESMF_FieldBundleAdd(sdat%pstrm(ns)%fldbun_model, (/lfield/), rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end do

       ! Create a field on the model mesh for coszen time interpolation for this stream if needed
       if (trim(sdat%stream(ns)%tinterpalgo) == 'coszen') then
          if (stream_nlev > 1) then
             sdat%pstrm(ns)%field_coszen = ESMF_FieldCreate(sdat%model_mesh, ESMF_TYPEKIND_r8, &
                  name='tavCosz', &
                  ungriddedLbound=(/1/), ungriddedUbound=(/stream_nlev/), gridToFieldMap=(/2/), &
                  meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
          else
             sdat%pstrm(ns)%field_coszen = ESMF_FieldCreate(sdat%model_mesh, ESMF_TYPEKIND_r8, &
                  name='tavCosz', &
                  meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
          end if
       endif

       ! ------------------------------------
       ! Create the sdat route handles for mapping the stream -> model
       ! ------------------------------------

       ! create the source and destination fields needed to create the route handles
       ! assume that all fields in a stream share the same mesh and there is only a unique model mesh
       ! can do this outside of a stream loop by just using the first stream index
       if(ESMF_MeshIsCreated(stream_mesh)) then
          if (stream_nlev > 1) then
             sdat%pstrm(ns)%field_stream = ESMF_FieldCreate(stream_mesh, &
                  ESMF_TYPEKIND_r8, meshloc=ESMF_MESHLOC_ELEMENT, &
                  ungriddedLbound=(/1/), ungriddedUbound=(/stream_nlev/), gridToFieldMap=(/2/), &
                  rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
          else
             sdat%pstrm(ns)%field_stream = ESMF_FieldCreate(stream_mesh, &
                  ESMF_TYPEKIND_r8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
             call ESMF_FieldFill(sdat%pstrm(ns)%field_stream, dataFillScheme="const", const1=1.0_r8, rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
          end if
       endif

       ! Why not use fldbun_model rather than fldbun_data?
       index = sdat%pstrm(ns)%stream_lb
       call dshr_fldbun_getFieldN(sdat%pstrm(ns)%fldbun_data(index), 1, lfield_dst, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       if (.not.  ESMF_MeshIsCreated(stream_mesh)) then
          sdat%stream(ns)%mapalgo = 'none'
       else
          if (trim(sdat%stream(ns)%mapalgo) == "bilinear") then
             call ESMF_FieldRegridStore(sdat%pstrm(ns)%field_stream, lfield_dst, &
                  routehandle=sdat%pstrm(ns)%routehandle, &
                  regridmethod=ESMF_REGRIDMETHOD_BILINEAR,  &
                  polemethod=ESMF_POLEMETHOD_ALLAVG, &
                  extrapMethod=ESMF_EXTRAPMETHOD_NEAREST_STOD, &
                  dstMaskValues = (/0/), &  ! ignore destination points where the mask is 0
                  srcMaskValues = (/0/), &  ! ignore source points where the mask is 0
                  srcTermProcessing=srcTermProcessing_Value, ignoreDegenerate=.true., rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
          else if (trim(sdat%stream(ns)%mapalgo) == 'redist') then
             call ESMF_FieldRedistStore(sdat%pstrm(ns)%field_stream, lfield_dst, &
                  routehandle=sdat%pstrm(ns)%routehandle, &
                  ignoreUnmatchedIndices = .true., rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
          else if (trim(sdat%stream(ns)%mapalgo) == 'nn') then
             call ESMF_FieldReGridStore(sdat%pstrm(ns)%field_stream, lfield_dst, &
                  routehandle=sdat%pstrm(ns)%routehandle, &
                  regridmethod=ESMF_REGRIDMETHOD_NEAREST_STOD, &
                  dstMaskValues = (/0/), &  ! ignore destination points where the mask is 0
                  srcMaskValues = (/0/), &  ! ignore source points where the mask is 0
                  srcTermProcessing=srcTermProcessing_Value, ignoreDegenerate=.true., &
                  unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)
          else if (trim(sdat%stream(ns)%mapalgo) == 'consf') then
             call ESMF_FieldReGridStore(sdat%pstrm(ns)%field_stream, lfield_dst, &
                  routehandle=sdat%pstrm(ns)%routehandle, &
                  regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
                  normType=ESMF_NORMTYPE_FRACAREA, &
                  dstMaskValues = (/0/), &  ! ignore destination points where the mask is 0
                  srcMaskValues = (/0/), &  ! ignore source points where the mask is 0
                  srcTermProcessing=srcTermProcessing_Value, ignoreDegenerate=.true., &
                  unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)
          else if (trim(sdat%stream(ns)%mapalgo) == 'consd') then
             call ESMF_FieldReGridStore(sdat%pstrm(ns)%field_stream, lfield_dst, &
                  routehandle=sdat%pstrm(ns)%routehandle, &
                  regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
                  normType=ESMF_NORMTYPE_DSTAREA, &
                  dstMaskValues = (/0/), &  ! ignore destination points where the mask is 0
                  srcMaskValues = (/0/), &  ! ignore source points where the mask is 0
                  srcTermProcessing=srcTermProcessing_Value, ignoreDegenerate=.true., &
                  unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)
          else if (trim(sdat%stream(ns)%mapalgo) == 'none') then
             ! single point stream data, no action required.
          else
             call shr_sys_abort('ERROR: map algo '//trim(sdat%stream(ns)%mapalgo)//' is not supported')
          end if
       end if

    end do ! end of loop over streams

    ! Check for vector pairs in the stream - BOTH ucomp and vcomp MUST BE IN THE SAME STREAM
    do ns = 1,shr_strdata_get_stream_count(sdat)
       stream_mesh => sdat%pstrm(ns)%stream_mesh
       stream_nlev = sdat%pstrm(ns)%stream_nlev
       stream_vector_names = trim(sdat%stream(ns)%stream_vectors)

       ! check that vector field list is a valid colon delimited string
       if (trim(stream_vector_names) /= 'null') then
          ! check that for now u and v are only for single level fields
          if (stream_nlev > 1) then
             ! TODO: add support for u and v for multi level fields
             call shr_sys_abort(subname//': vector fields are not currently supported for multi-level fields')
          end if
          ! check that stream vector names are valid
          if (.not. shr_string_listIsValid(stream_vector_names)) then
             write(sdat%stream(1)%logunit,*) trim(subname),' vec fldlist invalid m=',m,trim(stream_vector_names)
             call shr_sys_abort(subname//': vec fldlist invalid:'//trim(stream_vector_names))
          endif
          ! check that only 2 fields are contained for any vector pairing
          if (shr_string_listGetNum(stream_vector_names) /= 2) then
             write(sdat%stream(1)%logunit,*) trim(subname),' vec fldlist ne 2 m=',m,trim(stream_vector_names)
             call shr_sys_abort(subname//': vec fldlist ne 2:'//trim(stream_vector_names))
          endif
          ! create stream vector field
          sdat%pstrm(ns)%field_stream_vector = ESMF_FieldCreate(stream_mesh, &
               ESMF_TYPEKIND_r8, name='stream_vector', meshloc=ESMF_MESHLOC_ELEMENT, &
               ungriddedLbound=(/1/), ungriddedUbound=(/2/), gridToFieldMap=(/2/), rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
          if (mainproc) then
             write(sdat%stream(1)%logunit,'(a,i8)') "creating ESMF stream vector field with names" //&
                  trim(stream_vector_names)//" for stream ",ns
          end if
       end if
    enddo

    ! initialize sdat model clock and calendar
    sdat%model_clock = model_clock
    call ESMF_ClockGet(sdat%model_clock, calkindflag=esmf_caltype, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    if (esmf_caltype == ESMF_CALKIND_NOLEAP) then
       sdat%model_calendar = trim(shr_cal_noleap)
    else if (esmf_caltype == ESMF_CALKIND_GREGORIAN) then
       sdat%model_calendar = trim(shr_cal_gregorian)
    else
       call shr_sys_abort(subname//" ERROR bad ESMF calendar name "//trim(calendar))
    end if

    ! print sdat output
    if (mainproc) then
       if (present(stream_name)) then
          call shr_strdata_print(sdat, trim(stream_name))
       else
          call shr_strdata_print(sdat, 'stream_data')
       end if
       write(sdat%stream(1)%logunit,*) ' successfully initialized sdat'
    endif
  end subroutine shr_strdata_init

  !===============================================================================
  subroutine shr_strdata_get_stream_nlev(sdat, stream_index, rc)

    ! Obtain the number of vertical levels for the stream - this is obtained
    ! from the lev_dimname string

    ! input/output variables
    type(shr_strdata_type) , intent(inout) :: sdat
    integer                , intent(in)    :: stream_index
    integer                , intent(out)   :: rc

    ! local variables
    type(ESMF_VM)           :: vm
    type(file_desc_t)       :: pioid
    integer                 :: rcode
    character(CL)           :: filename
    integer                 :: dimid
    integer                 :: stream_nlev
    character(*), parameter :: subname = '(shr_strdata_set_stream_domain) '
    ! ----------------------------------------------

    rc = ESMF_SUCCESS

    ! Set ungridded dimension to the number of vertical level - must read this in
    if (trim(sdat%stream(stream_index)%lev_dimname) == 'null') then
       stream_nlev = 1
    else
       call ESMF_VMGetCurrent(vm, rc=rc)
       if (ChkErr(rc,__LINE__,u_FILE_u)) return
       if (sdat%mainproc) then
          call shr_stream_getData(sdat%stream(stream_index), 1, filename)
       end if
       call ESMF_VMBroadCast(vm, filename, CL, 0, rc=rc)
       rcode = pio_openfile(sdat%pio_subsystem, pioid, sdat%io_type, trim(filename), pio_nowrite)
       rcode = pio_inq_dimid(pioid, trim(sdat%stream(stream_index)%lev_dimname), dimid)
       rcode = pio_inq_dimlen(pioid, dimid, stream_nlev)
       call pio_closefile(pioid)
    end if
    if (sdat%mainproc) then
       write(sdat%stream(1)%logunit,*) trim(subname)//' stream_nlev = ',stream_nlev
    end if

    ! Set stream_nlev in the per-stream sdat info
    sdat%pstrm(stream_index)%stream_nlev = stream_nlev

  end subroutine shr_strdata_get_stream_nlev

  !===============================================================================
  subroutine shr_strdata_get_stream_domain(sdat, stream_index, fldname, flddata, rc)

    ! Obtain the data for fldname from the stream domain data

    ! input/output variables
    type(shr_strdata_type) , intent(inout) :: sdat
    integer                , intent(in)    :: stream_index
    character(len=*)       , intent(in)    :: fldname
    real(r8)               , pointer       :: flddata(:)
    integer                , intent(out)   :: rc

    ! local variables
    type(ESMF_VM)           :: vm
    type(var_desc_t)        :: varid
    type(file_desc_t)       :: pioid
    integer                 :: rcode
    character(CL)           :: filename
    type(io_desc_t)         :: pio_iodesc
    real(r4), allocatable   :: data_real(:)
    real(r8), allocatable   :: data_double(:)
    integer                 :: pio_iovartype
    integer                 :: lsize
    character(*), parameter :: subname = '(shr_strdata_set_stream_domain) '
    ! ----------------------------------------------

    rc = ESMF_SUCCESS
    call ESMF_VMGetCurrent(vm, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

    ! Determine the file to open
    if (sdat%mainproc) then
       call shr_stream_getData(sdat%stream(stream_index), 1, filename)
    end if
    call ESMF_VMBroadCast(vm, filename, CL, 0, rc=rc)

    ! Open the file
    rcode = pio_openfile(sdat%pio_subsystem, pioid, sdat%io_type, trim(filename), pio_nowrite)

    ! Create the pio iodesc for fldname
    call shr_strdata_set_stream_iodesc(sdat, sdat%pstrm(stream_index), trim(fldname), pioid, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

    ! Now read in the data for fldname
    lsize = size(flddata)
    rcode = pio_inq_varid(pioid, trim(fldname), varid)
    rcode = pio_inq_vartype(pioid, varid, pio_iovartype)
    if (pio_iovartype == PIO_REAL) then
       allocate(data_real(lsize))
       call pio_read_darray(pioid, varid, pio_iodesc, data_real, rcode)
       flddata(:) = real(data_real(:), kind=r8)
       deallocate(data_real)
    else if (pio_iovartype == PIO_DOUBLE) then
       allocate(data_double(lsize))
       call pio_read_darray(pioid, varid, pio_iodesc, data_double, rcode)
       flddata(:) = data_double(:)
       deallocate(data_double)
    else
       call shr_sys_abort(subName//"ERROR: only real and double types are supported for stream domain read")
    end if

    ! Free the memory associate with the iodesc and close the file
    call pio_freedecomp(pioid, pio_iodesc)
    call pio_closefile(pioid)

  end subroutine shr_strdata_get_stream_domain

  !===============================================================================
  subroutine shr_strdata_advance(sdat, ymd, tod, logunit, istr, timers, rc)

    ! -------------------------------------------------------
    ! Mismatching calendars: 4 cases
    ! (0) The stream calendar and model calendar are identical
    ! (1) The stream is a no leap calendar and the model is gregorian
    ! (2) The stream is a gregorian calendar and the model is a noleap calendar
    ! (3) The calendars mismatch and none of the above
    ! -------------------------------------------------------
    !
    ! ymdmod and todmod are the ymd and tod to time interpolate to.
    ! Generally, these are just the model date and time.  Also, always
    ! use the stream calendar for time interpolation for reasons
    ! described below.  When there is a calendar mismatch, support Feb
    ! 29 in a special way as needed to get reasonable values.  Note
    ! that when Feb 29 needs to be treated specially, a discontinuity
    ! will be introduced.  The size of that discontinuity will depend
    ! on the time series input data.
    !
    ! (0) The stream calendar and model calendar are identical:
    !     In this case it is still possible to have a mismatch if both are gregorian.
    !     These cases are:
    !     -   Model is no_leap, data is Gregorian and leapyear date 2/29 is encountered in data - skip date
    !     -   Model is Gregorian, data is no_leap and leapyear date 2/29 is encountered in model - repeat 2/28 data
    !     -   Model is Gregorian, data is gregorian but leapyears do not align.
    !     -       if in model leap year repeat data from 2/28
    !     -       if in data leap year skip date 2/29
    !
    !
    ! (1) The stream is a no leap calendar and the model is gregorian:
    ! Time interpolate on the noleap calendar.  If the model date is Feb 29,
    ! compute stream data for Feb 28 by setting ymdmod and todmod to Feb 28.
    ! This results in duplicate stream data on Feb 28 and Feb 29 and a
    ! discontinuity at the start of Feb 29.  This could potentially be fixed
    ! by using the gregorian calendar for time interpolation when the input data
    ! is relatively infrequent (say greater than daily) with the following concerns.
    !   - The forcing will not be reproduced identically on the same day with
    !     with climatological inputs data
    !   - Input data with variable input frequency might behave funny
    !   - An arbitrary discontinuity will be introduced in the time
    !     interpolation method based upon the logic chosen to transition
    !     from reproducing Feb 28 on Feb 29 and interpolating to Feb 29.
    !   - The time gradient of data will change by adding a day arbitrarily.
    !
    ! (2) The stream is a gregorian calendar and the model is a noleap calendar:
    ! Then just time interpolate on the gregorian calendar. This causes Feb 29
    ! stream data to be skipped and lead to a discontinuity at the start of March 1.
    !
    ! (3) If the calendars mismatch and neither of the three cases above
    ! are recognized, then abort.
    ! -------------------------------------------------------

    ! input/output variables
    type(shr_strdata_type) ,intent(inout)       :: sdat
    integer                ,intent(in)          :: ymd    ! current model date
    integer                ,intent(in)          :: tod    ! current model date
    integer                ,intent(in)          :: logunit
    character(len=*)       ,intent(in)          :: istr
    logical                ,intent(in),optional :: timers
    integer                ,intent(out)         :: rc

    ! local variables
    integer                             :: ns               ! stream index
    integer                             :: nf               ! field index
    integer                             :: i,lev,n          ! generic indices
    logical , allocatable               :: newData(:)
    integer , allocatable               :: ymdmod(:)        ! modified model dates to handle Feb 29
    real(r8), allocatable               :: coszen(:)        ! cosine of zenith angle
    integer                             :: todmod           ! modified model dates to handle Feb 29
    character(len=32)                   :: lstr             ! local string
    real(r8)                            :: flb,fub          ! factor for lb and ub
    real(r8) ,pointer                   :: dataptr1d(:)     ! pointer into field bundle
    real(r8) ,pointer                   :: dataptr1d_lb(:)  ! pointer into field bundle
    real(r8) ,pointer                   :: dataptr1d_ub(:)  ! pointer into field bundle
    real(r8) ,pointer                   :: dataptr2d(:,:)   ! pointer into field bundle
    real(r8) ,pointer                   :: dataptr2d_lb(:,:)! pointer into field bundle
    real(r8) ,pointer                   :: dataptr2d_ub(:,:)! pointer into field bundle
    real(r8), pointer                   :: nu_coords(:)     ! allocatable local element mesh lat and lons
    real(r8), pointer                   :: nv_coords(:)     ! allocatable local element mesh lat and lons
    real(r8), pointer                   :: data2d_src(:,:)  ! pointer into field bundle
    real(r8), pointer                   :: data2d_dst(:,:)  ! pointer into field bundle
    real(r8), pointer                   :: data_u_src(:)    ! pointer into field bundle
    real(r8), pointer                   :: data_v_src(:)    ! pointer into field bundle
    real(r8), pointer                   :: data_u_dst(:)    ! pointer into field bundle
    real(r8), pointer                   :: data_v_dst(:)    ! pointer into field bundle
    type(ESMF_Time)                     :: timeLB, timeUB   ! lb and ub times
    type(ESMF_TimeInterval)             :: timeint          ! delta time
    character(CL)                       :: calendar
    integer                             :: dday             ! delta days
    real(r8)                            :: dtime            ! delta time
    integer                             :: year,month,day   ! date year month day
    integer                             :: datayear,datamonth,dataday   ! data date year month day
    integer                             :: nstreams
    integer                             :: stream_index
    real(r8)         ,parameter         :: solZenMin = 0.001_r8 ! minimum solar zenith angle
    integer          ,parameter         :: tadj = 2
    character(len=*) ,parameter         :: timname = "_strd_adv"
    character(*)     ,parameter         :: subname = "(shr_strdata_advance) "
    character(*)     ,parameter         :: F00  = "('(shr_strdata_advance) ',a)"
    character(*)     ,parameter         :: F01  = "('(shr_strdata_advance) ',a,a,i4,2(f10.5,2x))"
    !-------------------------------------------------------------------------------

    rc = ESMF_SUCCESS

    nullify(dataptr1d)
    nullify(dataptr1d_ub)
    nullify(dataptr1d_lb)
    nullify(dataptr2d)
    nullify(dataptr2d_ub)
    nullify(dataptr2d_lb)
    nullify(nu_coords)
    nullify(nv_coords)
    nullify(data2d_src)
    nullify(data2d_dst)
    nullify(data_u_src)
    nullify(data_v_src)
    nullify(data_u_dst)
    nullify(data_v_dst)

    nstreams = shr_strdata_get_stream_count(sdat)
    if (nstreams < 1) return ! TODO: is this needed

    lstr = trim(istr)
    ! To avoid an unused dummy variable warning
    if(present(timers)) then
       write(sdat%stream(1)%logunit,*) trim(subname),'optional variable timers present but unused'
    endif
    call ESMF_TraceRegionEnter(trim(lstr)//trim(timname)//'_total')

    sdat%ymd = ymd
    sdat%tod = tod
    if (nstreams > 0) then
       allocate(newData(nstreams))
       allocate(ymdmod(nstreams))

       do ns = 1,nstreams
          ! ---------------------------------------------------------
          ! Consistency checks
          ! ---------------------------------------------------------

          ! case(0)
          ymdmod(ns) = ymd
          todmod    = tod
          calendar = trim(sdat%stream(ns)%calendar)
          call shr_cal_date2ymd (ymd,year,month,day)
          if (trim(sdat%model_calendar) /= trim(sdat%stream(ns)%calendar)) then
             if (( trim(sdat%model_calendar) == trim(shr_cal_gregorian)) .and. &
                  (trim(sdat%stream(ns)%calendar) == trim(shr_cal_noleap))) then
                ! case (1), set feb 29 = feb 28
                call shr_cal_date2ymd (ymd,year,month,day)
                if (month == 2 .and. day == 29) then
                   call shr_cal_ymd2date(year,2,28,ymdmod(ns))
                endif
                calendar = shr_cal_noleap
             endif
          endif
          ! ---------------------------------------------------------
          ! Determine if new data is read in - if so then copy
          ! fldbun_stream_ub to fldbun_stream_lb and read in new fldbun_stream_ub data
          ! ---------------------------------------------------------

          call ESMF_TraceRegionEnter(trim(lstr)//trim(timname)//'_readLBUB')

          select case(sdat%stream(ns)%readmode)
          case ('single')
             call shr_strdata_readLBUB(sdat, ns, ymdmod(ns), todmod, newData(ns), trim(lstr)//'_readLBUB', rc=rc)
             if (ChkErr(rc,__LINE__,u_FILE_u)) return
          case ('full_file')
             ! TODO: need to put in capability to read all stream data at once
          case default
             write(logunit,F00) "ERROR: Unsupported readmode : ", trim(sdat%stream(ns)%readmode)
             call shr_sys_abort(subName//"ERROR: Unsupported readmode: "//trim(sdat%stream(ns)%readmode))
          end select

          if (debug > 0 .and. sdat%mainproc) then
             write(sdat%stream(1)%logunit,*) trim(subname),' newData flag = ',ns,newData(ns)
             write(sdat%stream(1)%logunit,*) trim(subname),' LB ymd,tod = ',ns,sdat%pstrm(ns)%ymdLB,sdat%pstrm(ns)%todLB
             write(sdat%stream(1)%logunit,*) trim(subname),' UB ymd,tod = ',ns,sdat%pstrm(ns)%ymdUB,sdat%pstrm(ns)%todUB
          endif

          ! ---------------------------------------------------------
          ! If new data is read in:
          ! ---------------------------------------------------------

          if (newData(ns)) then
             if (trim(sdat%model_calendar) /= trim(sdat%stream(ns)%calendar)) then
                if ((trim(sdat%model_calendar) == trim(shr_cal_noleap)) .and. &
                     (trim(sdat%stream(ns)%calendar) == trim(shr_cal_gregorian))) then
                   call shr_cal_date2ymd(sdat%pstrm(ns)%ymdUB, datayear, datamonth, dataday)
                   if(datamonth==3 .and. dataday==1 .and. month==2 .and. day==28) then
                      calendar = shr_cal_noleap
                   endif
                   ! case (2), feb 29 input data will be skipped automatically
                else if (.not. ( trim(sdat%model_calendar) == trim(shr_cal_gregorian)) .and. &
                  (trim(sdat%stream(ns)%calendar) == trim(shr_cal_noleap))) then
                   ! case (3), abort
                   write(logunit,*) trim(subname),' ERROR: mismatch calendar ', &
                        trim(sdat%model_calendar),':',trim(sdat%stream(ns)%calendar)
                   call shr_sys_abort(trim(subname)//' ERROR: mismatch calendar ')
                endif
             else ! calendars are the same
                if(trim(sdat%model_calendar) == trim(shr_cal_gregorian)) then
                   ! Both are in gregorian - but it's possible that there is a mismatch
                   ! such that the model is in leapyear but the data is not
                   call shr_cal_date2ymd (ymd,year,month,day)
                   call shr_cal_date2ymd(sdat%pstrm(ns)%ymdUB, datayear, datamonth, dataday)
                   if(month == 2 .and. day >= 28) then
                      if(shr_cal_leapyear(year) .and. .not. shr_cal_leapyear(datayear)) then
                         ! model is in leap year but data is not
                         calendar = shr_cal_noleap
                      endif
                   else
                      calendar = sdat%model_calendar
                   endif
                else
                   calendar = sdat%model_calendar
                endif
             endif

             ! Reset time bounds if newdata read in
             call shr_cal_timeSet(timeLB,sdat%pstrm(ns)%ymdLB,0,calendar,rc=rc)
             if (ChkErr(rc,__LINE__,u_FILE_u)) return
             call shr_cal_timeSet(timeUB,sdat%pstrm(ns)%ymdUB,0,calendar,rc=rc)
             if (ChkErr(rc,__LINE__,u_FILE_u)) return

             timeint = timeUB-timeLB
             call ESMF_TimeIntervalGet(timeint, StartTimeIn=timeLB, d=dday)
             if (ChkErr(rc,__LINE__,u_FILE_u)) return

             if (.not. sdat%pstrm(ns)%override_annual_cycle) then
                if(sdat%stream(ns)%dtlimit == -1) then
                   sdat%pstrm(ns)%override_annual_cycle = .true.
                   if(sdat%mainproc) then
                      write(logunit,*) trim(subname),' WARNING: Stream ',ns,' is not cycling on annual boundaries, and dtlimit check has been overridden'
                   endif
                else
                   dtime = abs(real(dday,r8) + real(sdat%pstrm(ns)%todUB-sdat%pstrm(ns)%todLB,r8)/shr_const_cDay)

                   sdat%pstrm(ns)%dtmin = min(sdat%pstrm(ns)%dtmin,dtime)
                   sdat%pstrm(ns)%dtmax = max(sdat%pstrm(ns)%dtmax,dtime)

                   if ((sdat%pstrm(ns)%dtmax/sdat%pstrm(ns)%dtmin) > sdat%stream(ns)%dtlimit) then
                      if (sdat%mainproc) then
                         write(sdat%stream(1)%logunit,'(a,i8)') trim(subname),' ERROR: for stream ',ns
                         write(sdat%stream(1)%logunit,'(a,i8)') trim(subname),' ERROR: dday = ',dday
                         write(sdat%stream(1)%logunit,'(a,4(f15.5,2x))') trim(subName),' ERROR: dtime, dtmax, dtmin, dtlimit = ',&
                              dtime, sdat%pstrm(ns)%dtmax, sdat%pstrm(ns)%dtmin, sdat%stream(ns)%dtlimit
                         write(sdat%stream(1)%logunit,'(a,4(i10,2x))') trim(subName),' ERROR: ymdLB, todLB, ymdUB, todUB = ', &
                              sdat%pstrm(ns)%ymdLB, sdat%pstrm(ns)%todLB, sdat%pstrm(ns)%ymdUB, sdat%pstrm(ns)%todUB
                      end if
                      write(6,*) trim(subname),' ERROR: for stream ',ns, ' and calendar ',trim(calendar)
                      write(6,*) trim(subName),' ERROR: dtime, dtmax, dtmin, dtlimit = ',&
                           dtime, sdat%pstrm(ns)%dtmax, sdat%pstrm(ns)%dtmin, sdat%stream(ns)%dtlimit
                      write(6,*) trim(subName),' ERROR: ymdLB, todLB, ymdUB, todUB = ', &
                           sdat%pstrm(ns)%ymdLB, sdat%pstrm(ns)%todLB, sdat%pstrm(ns)%ymdUB, sdat%pstrm(ns)%todUB
                      call shr_sys_abort(trim(subName)//' ERROR dt limit for stream, see atm.log output')
                   endif
                endif
             endif
          endif

          call ESMF_TraceRegionExit(trim(lstr)//trim(timname)//'_readLBUB')

       enddo ! end of loop over streams

       ! ---------------------------------------------------------
       ! Do time interpolation to create fldbun_model
       ! ---------------------------------------------------------

       do ns = 1,nstreams
          if (trim(sdat%stream(ns)%tinterpalgo) == 'coszen') then

             ! Determine stream lower bound index
             stream_index = sdat%pstrm(ns)%stream_lb

             ! ------------------------------------------
             ! time interpolation method is coszen
             ! ------------------------------------------

             call ESMF_TraceRegionEnter(trim(lstr)//trim(timname)//'_coszen')
             allocate(coszen(sdat%model_lsize))

             ! get coszen
             call ESMF_TraceRegionEnter(trim(lstr)//trim(timname)//'_coszenC')
             call shr_tInterp_getCosz(coszen, sdat%model_lon, sdat%model_lat, ymdmod(ns), todmod, &
                  sdat%eccen, sdat%mvelpp, sdat%lambm0, sdat%obliqr, sdat%stream(ns)%calendar, &
                  sdat%mainproc, sdat%stream(1)%logunit)
             call ESMF_TraceRegionExit(trim(lstr)//trim(timname)//'_coszenC')
             if (debug > 0 .and. sdat%mainproc) then
                do n = 1,size(coszen)
                   write(sdat%stream(1)%logunit,'(a,i4,2x,2(i18,2x),i8,d20.10)')' stream,ymdmod,todmod,n,coszen= ',&
                        ns, ymd, tod, n, coszen(n)
                end do
             end if

             ! get avg cosz factor
             if (newdata(ns)) then
                ! compute a new avg cosz
                call ESMF_TraceRegionEnter(trim(lstr)//trim(timname)//'_coszenN')
                if (.not. allocated(sdat%tavCoszen)) then
                   allocate(sdat%tavCoszen(sdat%model_lsize))
                end if
                call shr_tInterp_getAvgCosz(sdat%tavCoszen, sdat%model_lon, sdat%model_lat,  &
                     sdat%pstrm(ns)%ymdLB, sdat%pstrm(ns)%todLB,  sdat%pstrm(ns)%ymdUB, sdat%pstrm(ns)%todUB,  &
                     sdat%eccen, sdat%mvelpp, sdat%lambm0, sdat%obliqr, sdat%modeldt, &
                     sdat%stream(ns)%calendar, sdat%mainproc, sdat%stream(1)%logunit, rc=rc)
                call ESMF_TraceRegionExit(trim(lstr)//trim(timname)//'_coszenN')
                if (debug > 0 .and. sdat%mainproc) then
                   do n = 1,size(coszen)
                      write(sdat%stream(1)%logunit,'(a,i4,2x,4(i18,2x),i8,d20.10)')' stream,lbymd,lbsec,ubymd,ubsec,newdata,n,tavgCoszen= ',&
                           ns, sdat%pstrm(ns)%ymdLB, sdat%pstrm(ns)%todLB, sdat%pstrm(ns)%ymdUB, sdat%pstrm(ns)%todUB, &
                           n, sdat%tavCoszen(n)
                   end do
                end if
             endif

             ! compute time interperpolate value - LB data normalized with this factor: cosz/tavCosz
             do nf = 1,size(sdat%pstrm(ns)%fldlist_model)
                if (sdat%pstrm(ns)%stream_nlev > 1) then
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_model, &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr2=dataptr2d, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_data(stream_index), &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr2=dataptr2d_lb, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   do i = 1,size(dataptr2d,dim=2)
                      if (coszen(i) > solZenMin) then
                         dataptr2d(:,i) = dataptr2d_lb(:,i)*coszen(i)/sdat%tavCoszen(i)
                      else
                         dataptr2d(:,i) = 0._r8
                      endif
                   end do
                else
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_model, &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr1=dataptr1d, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_data(stream_index), &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr1=dataptr1d_lb, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   do i = 1,size(dataptr1d)
                      if (coszen(i) > solZenMin) then
                         dataptr1d(i) = dataptr1d_lb(i)*coszen(i)/sdat%tavCoszen(i)
                      else
                         dataptr1d(i) = 0._r8
                      endif
                   end do
                end if
             end do

             deallocate(coszen)
             call ESMF_TraceRegionExit(trim(lstr)//trim(timname)//'_coszen')

          elseif (trim(sdat%stream(ns)%tinterpalgo) /= trim(shr_strdata_nullstr)) then

             ! ------------------------------------------
             ! time interpolation method is not coszen
             ! ------------------------------------------

             call ESMF_TraceRegionEnter(trim(lstr)//trim(timname)//'_tint')
             call shr_tInterp_getFactors(sdat%pstrm(ns)%ymdlb, sdat%pstrm(ns)%todlb, &
                  sdat%pstrm(ns)%ymdub, sdat%pstrm(ns)%todub, &
                  ymdmod(ns), todmod, flb, fub, calendar=sdat%stream(ns)%calendar, logunit=sdat%stream(1)%logunit, &
                  algo=trim(sdat%stream(ns)%tinterpalgo), rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
             if (debug > 0 .and. sdat%mainproc) then
                write(sdat%stream(1)%logunit,'(a,i4,2(f10.5,2x))') &
                     trim(subname)//' non-cosz-interp stream, flb, fub= ',ns,flb,fub
             endif
             do nf = 1,size(sdat%pstrm(ns)%fldlist_model)
                if (sdat%pstrm(ns)%stream_nlev > 1) then
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_model, &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr2=dataptr2d, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_data(sdat%pstrm(ns)%stream_lb), &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr2=dataptr2d_lb, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_data(sdat%pstrm(ns)%stream_ub), &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr2=dataptr2d_ub, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   do lev = 1,sdat%pstrm(ns)%stream_nlev
                      dataptr2d(lev,:) = dataptr2d_lb(lev,:) * flb + dataptr2d_ub(lev,:) * fub
                   end do
                else
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_model, &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr1=dataptr1d, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_data(sdat%pstrm(ns)%stream_lb), &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr1=dataptr1d_lb, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_data(sdat%pstrm(ns)%stream_ub), &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr1=dataptr1d_ub, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   dataptr1d(:) = dataptr1d_lb(:) * flb + dataptr1d_ub(:) * fub
                end if
             end do
             call ESMF_TraceRegionExit(trim(lstr)//trim(timname)//'_tint')

          else

             ! ------------------------------------------
             ! zero out stream data for this field
             ! ------------------------------------------

             call ESMF_TraceRegionEnter(trim(lstr)//trim(timname)//'_zero')
             do nf = 1,size(sdat%pstrm(ns)%fldlist_model)
                if (sdat%pstrm(ns)%stream_nlev > 1) then
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_model, &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr2=dataptr2d, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   dataptr2d(:,:) = 0._r8
                else
                   call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_model, &
                        sdat%pstrm(ns)%fldlist_model(nf), fldptr1=dataptr1d, rc=rc)
                   if (ChkErr(rc,__LINE__,u_FILE_u)) return
                   dataptr1d(:) = 0._r8
                end if
             end do
             call ESMF_TraceRegionExit(trim(lstr)//trim(timname)//'_zero')

          endif
       end do  ! loop over ns (number of streams)

       deallocate(newData)
       deallocate(ymdmod)

    endif  ! nstreams > 0

    call ESMF_TraceRegionExit(trim(lstr)//trim(timname)//'_total')

  end subroutine shr_strdata_advance

  !===============================================================================

  subroutine shr_strdata_setOrbs(sdat,eccen,mvelpp,lambm0,obliqr,modeldt)

    type(shr_strdata_type),intent(inout) :: sdat
    real(r8),intent(in) :: eccen
    real(r8),intent(in) :: mvelpp
    real(r8),intent(in) :: lambm0
    real(r8),intent(in) :: obliqr
    integer,intent(in) :: modeldt

    ! local variables
    character(len=*),parameter :: subname = "(shr_strdata_setOrbs) "
    !-------------------------------------------------------------------------------

    sdat%eccen   = eccen
    sdat%mvelpp  = mvelpp
    sdat%lambm0  = lambm0
    sdat%obliqr  = obliqr
    sdat%modeldt = modeldt

  end subroutine shr_strdata_setOrbs

  !===============================================================================
  subroutine shr_strdata_print(sdat, name)

    !  Print strdata common to all data models

    ! input/output parameters:
    type(shr_strdata_type) , intent(in) :: sdat  ! strdata data data-type
    character(len=*)       , intent(in) :: name

    ! local variables
    integer  :: ns
    character(*),parameter :: subName = "(shr_strdata_print) "
    character(*),parameter ::   F00 = "('(shr_strdata_print) ',8a)"
    character(*),parameter ::   F01 = "('(shr_strdata_print) ',a,i6,a)"
    character(*),parameter ::   F02 = "('(shr_strdata_print) ',a,es13.6)"
    character(*),parameter ::   F04 = "('(shr_strdata_print) ',a,i2,a,a)"
    character(*),parameter ::   F05 = "('(shr_strdata_print) ',a)"
    character(*),parameter ::   F07 = "('(shr_strdata_print) ',a,i2,a,es13.6)"
    character(*),parameter ::   F90 = "('(shr_strdata_print) ',58('-'))"
    !-------------------------------------------------------------------------------

    write(sdat%stream(1)%logunit,*)
    write(sdat%stream(1)%logunit,F90)
    write(sdat%stream(1)%logunit,F00) "name        = ",trim(name)
    write(sdat%stream(1)%logunit,F00) "calendar    = ",trim(sdat%model_calendar)
    write(sdat%stream(1)%logunit,F02) "eccen       = ",sdat%eccen
    write(sdat%stream(1)%logunit,F02) "mvelpp      = ",sdat%mvelpp
    write(sdat%stream(1)%logunit,F02) "lambm0      = ",sdat%lambm0
    write(sdat%stream(1)%logunit,F02) "obliqr      = ",sdat%obliqr
    write(sdat%stream(1)%logunit,F01) "pio_iotype  = ",sdat%io_type
    write(sdat%stream(1)%logunit,F01) "nstreams    = ",shr_strdata_get_stream_count(sdat)
    write(sdat%stream(1)%logunit,F05) "Per stream information "
    do ns = 1, shr_strdata_get_stream_count(sdat)
       write(sdat%stream(1)%logunit,F04) "  taxMode (",ns,") = ",trim(sdat%stream(ns)%taxmode)
       write(sdat%stream(1)%logunit,F07) "  dtlimit (",ns,") = ",sdat%stream(ns)%dtlimit
       write(sdat%stream(1)%logunit,F04) "  mapalgo (",ns,") = ",trim(sdat%stream(ns)%mapalgo)
       write(sdat%stream(1)%logunit,F04) "  tintalgo(",ns,") = ",trim(sdat%stream(ns)%tinterpalgo)
       write(sdat%stream(1)%logunit,F04) "  readmode(",ns,") = ",trim(sdat%stream(ns)%readmode)
       write(sdat%stream(1)%logunit,F04) "  vectors (",ns,") = ",trim(sdat%stream(ns)%stream_vectors)
       write(sdat%stream(1)%logunit,F01) " "
    end do
    write(sdat%stream(1)%logunit,F90)

  end subroutine shr_strdata_print

  !===============================================================================

  subroutine shr_strdata_readLBUB(sdat, ns, mDate, mSec, newData, istr, rc)

    !-------------------------------------------------------------------------
    ! Read LB and UB of stream data
    !-------------------------------------------------------------------------

    ! input/output variables
    type(shr_strdata_type) , intent(inout), target :: sdat  ! strdata data data-type
    integer                , intent(in)            :: ns    ! stream number
    integer                , intent(in)            :: mDate, mSec
    logical                , intent(out)           :: newData
    character(len=*)       , intent(in)            :: istr
    integer                , intent(out)           :: rc

    ! local variables
    type(shr_stream_streamType), pointer :: stream
    type(ESMF_Mesh)            , pointer :: stream_mesh
    type(ESMF_VM)                        :: vm
    logical                              :: fileexists
    integer                              :: oDateLB,oSecLB,dDateLB
    integer                              :: oDateUB,oSecUB,dDateUB
    real(r8)                             :: rDateM,rDateLB,rDateUB  ! model,LB,UB dates with fractional days
    integer                              :: n_lb, n_ub
    integer                              :: i
    character(CL)                        :: filename_lb
    character(CL)                        :: filename_ub
    character(CL)                        :: filename_next
    character(CL)                        :: filename_prev
    logical                              :: find_bounds
    character(*), parameter              :: subname = '(shr_strdata_readLBUB) '
    character(*), parameter              :: F00   = "('(shr_strdata_readLBUB) ',8a)"
    character(*), parameter              :: F01   = "('(shr_strdata_readLBUB) ',a,5i8)"
    !-------------------------------------------------------------------------

    rc = ESMF_SUCCESS
    stream => sdat%stream(ns)
    stream_mesh => sdat%pstrm(ns)%stream_mesh

    call ESMF_TraceRegionEnter(trim(istr)//'_setup')
    ! allocate streamdat instance on all tasks
    call ESMF_VMGetCurrent(vm, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return

    newData = .false.
    n_lb = -1
    n_ub = -1
    filename_lb = 'undefinedlb'
    filename_ub = 'undefinedub'

    ! stream lower bound and upper bound dates
    oDateLB = sdat%pstrm(ns)%ymdLB
    oSecLB  = sdat%pstrm(ns)%todLB
    oDateUB = sdat%pstrm(ns)%ymdUB
    oSecUB  = sdat%pstrm(ns)%todUB

    ! module current, lowerbound and upperbound date
    rDateM  = real(mDate  ,r8) + real(mSec  ,r8)/shr_const_cday
    rDateLB = real(sdat%pstrm(ns)%ymdLB,r8) + real(sdat%pstrm(ns)%todLB,r8)/shr_const_cday
    rDateUB = real(sdat%pstrm(ns)%ymdUB,r8) + real(sdat%pstrm(ns)%todUB,r8)/shr_const_cday

    call ESMF_TraceRegionExit(trim(istr)//'_setup')

    ! if model current date is outside of model lower or upper bound - find the stream bounds
    find_bounds = (rDateM < rDateLB .or. rDateM >= rDateUB)
    if (debug > 0 .and. sdat%mainproc) then
       write(sdat%stream(1)%logunit,'(a,i4,2x,6(i18,2x),l7)')' stream,lbymd,lbsec,mdate,msec,ubymd,ubsec,newdata= ',ns,&
            sdat%pstrm(ns)%ymdLB,sdat%pstrm(ns)%todLB, &
            mdate,msec, &
            sdat%pstrm(ns)%ymdUB,sdat%pstrm(ns)%todUB,find_bounds
       write(sdat%stream(1)%logunit,'(a,i4,2x,3(f20.3,2x),l7)')' stream,rdateLB,rdateM,rdateUB,newdata= ',&
            ns,rdateLB,rdateM,rdateUB,find_bounds
    end if
    if (find_bounds) then
       call ESMF_TraceRegionEnter(trim(istr)//'_fbound')
       call shr_stream_findBounds(stream, mDate, mSec,  sdat%mainproc, &
            sdat%pstrm(ns)%ymdLB, dDateLB, sdat%pstrm(ns)%todLB, n_lb, filename_lb, &
            sdat%pstrm(ns)%ymdUB, dDateUB, sdat%pstrm(ns)%todUB, n_ub, filename_ub)
       call ESMF_TraceRegionExit(trim(istr)//'_fbound')
       if (debug > 0 .and. sdat%mainproc) then
          write(sdat%stream(1)%logunit,'(a,i4,2x,6(i18,2x),l7)')' stream,lbymd,lbsec,mdate,msec,ubymd,ubsec,newdata= ',ns,&
               sdat%pstrm(ns)%ymdLB,sdat%pstrm(ns)%todLB,&
               mdate,msec, &
               sdat%pstrm(ns)%ymdUB,sdat%pstrm(ns)%todUB
          write(sdat%stream(1)%logunit,'(a,i4,2x,3(f20.3,2x),l7)')' stream,rdateLB,rdateM,rdateUB,newdata= ',&
               ns,rdateLB,rdateM,rdateUB,find_bounds
       end if
    endif

    ! determine if need to read in new stream data
    newdata = (sdat%pstrm(ns)%ymdLB /= oDateLB .or. sdat%pstrm(ns)%todLB /= oSecLB)
    if (newdata) then
       if (sdat%pstrm(ns)%ymdLB == oDateUB .and. sdat%pstrm(ns)%todLB == oSecUB) then
          ! copy fldbun_stream_ub to fldbun_stream_lb
          i = sdat%pstrm(ns)%stream_ub
          sdat%pstrm(ns)%stream_ub = sdat%pstrm(ns)%stream_lb
          sdat%pstrm(ns)%stream_lb = i
       else
          ! read lower bound of data
          call shr_strdata_readstrm(sdat, sdat%pstrm(ns), stream, &
               sdat%pstrm(ns)%fldbun_data(sdat%pstrm(ns)%stream_lb), &
               filename_lb, n_lb, istr=trim(istr)//'_LB', boundstr='lb', rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if
    endif

    if (newdata) then
       ! read upper bound of data
       call shr_strdata_readstrm(sdat, sdat%pstrm(ns), stream, &
            sdat%pstrm(ns)%fldbun_data(sdat%pstrm(ns)%stream_ub), &
            filename_ub, n_ub, istr=trim(istr)//'_UB', boundstr='ub', rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
    endif

    ! determine previous & next data files in list of files
    call ESMF_TraceRegionEnter(trim(istr)//'_filemgt')
    if (sdat%mainproc .and. newdata) then
       call shr_stream_getPrevFileName(stream, filename_lb, filename_prev)
       call shr_stream_getNextFileName(stream, filename_ub, filename_next)
       inquire(file=trim(filename_next),exist=fileExists)
       if ( trim(filename_next) == "unknown" .or. fileExists) then
          ! do nothing
       end if
    endif
    call ESMF_TraceRegionExit(trim(istr)//'_filemgt')

  end subroutine shr_strdata_readLBUB

  !===============================================================================
  subroutine shr_strdata_readstrm(sdat, per_stream, stream, fldbun_data, &
       filename, nt, istr, boundstr, rc)

    ! Read the stream data and initialize the strea pio_iodesc the first time
    ! the stream is read

    use shr_const_mod         , only : r8fill => SHR_CONST_SPVAL
    use shr_infnan_mod        , only : shr_infnan_isnan

    ! input/output variables
    type(shr_strdata_type)      , intent(inout)         :: sdat  ! strdata data data-type
    type(shr_strdata_perstream) , intent(inout)         :: per_stream
    type(shr_stream_streamType) , intent(inout)         :: stream
    type(ESMF_FieldBundle)      , intent(inout)         :: fldbun_data
    character(len=*)            , intent(in)            :: filename
    integer                     , intent(in)            :: nt
    character(len=*)            , intent(in)            :: istr
    character(len=*)            , intent(in)            :: boundstr
    integer                     , intent(out)           :: rc

    ! local variables
    integer                  :: stream_nlev
    type(ESMF_Field)         :: field_dst, field_vector_dst
    character(CL)            :: currfile
    logical                  :: fileexists
    logical                  :: fileopen
    type(file_desc_t)        :: pioid
    type(var_desc_t)         :: varid
    integer                  :: nf
    integer                  :: rCode
    real(r4)                 :: fillvalue_r4
    real(r8)                 :: fillvalue_r8
    logical                  :: handlefill = .false.
    integer                  :: old_error_handle
    real(r8), pointer        :: dataptr(:)
    real(r8), pointer        :: dataptr1d(:)        ! field bundle data
    real(r8), pointer        :: dataptr2d(:,:)      ! field bundle data
    real(r8), pointer        :: dataptr2d_src(:,:)  ! field bundle data
    real(r8), pointer        :: dataptr2d_dst(:,:)  ! field bundle data
    real(r4), allocatable    :: data_real1d(:)      ! stream input data
    real(r4), allocatable    :: data_real2d(:,:)    ! stream input data
    real(r8), allocatable    :: data_dbl1d(:)       ! stream input data
    real(r8), allocatable    :: data_dbl2d(:,:)     ! stream input data
    integer(i2), allocatable :: data_short1d(:)     ! stream input data
    integer(i2), allocatable :: data_short2d(:,:)   ! stream input data
    integer                  :: lsize, n
    integer                  :: spatialDim, numOwnedElements
    integer                  :: pio_iovartype
    real(r8), pointer        :: nv_coords(:), nu_coords(:)
    real(r8), pointer        :: data_u_dst(:), data_v_dst(:)
    real(r8)                 :: lat, lon
    real(r8)                 :: sinlat, sinlon
    real(r8)                 :: coslat, coslon
    real(r8)                 :: scale_factor, add_offset
    integer(i2)              :: fillvalue_i2
    character(CS)            :: uname, vname
    integer                  :: i, lev
    logical                  :: checkflag = .false.
    character(*), parameter  :: subname = '(shr_strdata_readstrm) '
    character(*), parameter  :: F00   = "('(shr_strdata_readstrm) ',8a)"
    character(*), parameter  :: F02   = "('(shr_strdata_readstrm) ',2a,i8)"
    character(CL)            :: errmsg
    !-------------------------------------------------------------------------------

    rc = ESMF_SUCCESS

    ! nullify local pointers
    nullify(dataptr)
    nullify(dataptr1d)
    nullify(dataptr2d)
    nullify(dataptr2d_src)
    nullify(dataptr2d_dst)
    nullify(nv_coords)
    nullify(nu_coords)
    nullify(data_u_dst)
    nullify(data_v_dst)

    ! Set up file to read from
    if (sdat%mainproc) then
       inquire(file=trim(fileName),exist=fileExists)
       if (.not. fileExists) then
          write(sdat%stream(1)%logunit,F00) "ERROR: file does not exist: ", trim(fileName)
          call shr_sys_abort(subName//"ERROR: file does not exist: "//trim(fileName))
       end if
    endif

    ! Get current file and determine if it is open
    call shr_stream_getCurrFile(stream, fileopen=fileopen, currfile=currfile, currpioid=pioid)

    if (fileopen .and. currfile==filename) then
       ! don't reopen file, all good
    else
       ! otherwise close the old file if open and open new file
       if (fileopen) then
          if (sdat%mainproc) write(sdat%stream(1)%logunit,F00) 'close  : ',trim(currfile)
          call pio_closefile(pioid)
       endif
       if (sdat%mainproc) write(sdat%stream(1)%logunit,F00) 'opening   : ',trim(filename)
       rcode = pio_openfile(sdat%pio_subsystem, pioid, sdat%io_type, trim(filename), pio_nowrite)
       call shr_stream_setCurrFile(stream, fileopen=.true., currfile=trim(filename), currpioid=pioid)
    endif

    ! ******************************************************************************
    ! Determine the pio io descriptor for the stream from the first data field in the stream
    ! ******************************************************************************

    stream_nlev = per_stream%stream_nlev

    if (ESMF_MeshIsCreated(per_stream%stream_mesh)) then
       if (.not. per_stream%stream_pio_iodesc_set) then
          if (sdat%mainproc) write(sdat%stream(1)%logunit,F00) 'setting pio descriptor : ',trim(filename)
          call shr_strdata_set_stream_iodesc(sdat, per_stream, trim(per_stream%fldlist_stream(1)), &
               pioid, rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
          per_stream%stream_pio_iodesc_set = .true.
       endif
       if (stream_nlev > 1) then
          call dshr_field_getfldptr(per_stream%field_stream, fldptr2=dataptr2d, rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       else
          call dshr_field_getfldptr(per_stream%field_stream, fldptr1=dataptr1d, rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if
    else
       if (stream_nlev == 1) then
          allocate(dataptr1d(1))
       else
          call shr_sys_abort("ERROR: multi-level streams always require a stream mesh")
       end if
    end if

    ! ******************************************************************************
    ! Read in the stream data for field names in fldname_stream_input - but fill in
    ! the data for fldbun_stream with the field names fldname_stream_model
    ! Note that the dimensions in the field bundle field are (stream_nlev,i) whereas
    ! the dimension of the input data are (i,stream_nlev)
    ! ******************************************************************************

    call ESMF_TraceRegionEnter(trim(istr)//'_readpio')
    if (sdat%mainproc) then
       write(sdat%stream(1)%logunit,F02) 'reading file ' // trim(boundstr) //': ',trim(filename), nt
    endif

    if (ESMF_FieldIsCreated(per_stream%field_stream_vector)) then
       call shr_string_listGetName(stream%stream_vectors,1,uname)
       call shr_string_listGetName(stream%stream_vectors,2,vname)
       call dshr_field_getfldptr(per_stream%field_stream_vector, fldptr2=dataptr2d_src, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
    endif

    do nf = 1,size(per_stream%fldlist_stream)

       ! determine type of the variable
       rcode = pio_inq_varid(pioid, trim(per_stream%fldlist_stream(nf)), varid)
       rcode = pio_inq_vartype(pioid, varid, pio_iovartype)

       ! allocate memory for input read
       if (stream_nlev > 1) then
          lsize = size(dataptr2d, dim=2)
          if (pio_iovartype == PIO_REAL .and. .not. allocated(data_real2d)) then
             allocate(data_real2d(lsize, stream_nlev))
          else if (pio_iovartype == PIO_DOUBLE .and. .not. allocated(data_dbl2d)) then
             allocate(data_dbl2d(lsize, stream_nlev))
          else if(pio_iovartype == PIO_SHORT .and. .not. allocated(data_short2d)) then
             allocate(data_short2d(lsize, stream_nlev))
          endif
       else
          lsize = size(dataptr1d)
          if (pio_iovartype == PIO_REAL .and. .not. allocated(data_real1d)) then
             allocate(data_real1d(lsize))
          else if (pio_iovartype == PIO_DOUBLE .and. .not. allocated(data_dbl1d)) then
             allocate(data_dbl1d(lsize))
          else if(pio_iovartype == PIO_SHORT .and. .not. allocated(data_short1d)) then
             allocate(data_short1d(lsize))
          endif
       end if

       ! determine if will handle fill
       handlefill = .false.
       call PIO_seterrorhandling(pioid, PIO_BCAST_ERROR, old_error_handle)
       if (pio_iovartype == PIO_REAL) then
          rcode = pio_get_att(pioid, varid, "_FillValue", fillvalue_r4)
       else if (pio_iovartype == PIO_DOUBLE) then
          rcode = pio_get_att(pioid, varid, "_FillValue", fillvalue_r8)
       else if (pio_iovartype == PIO_SHORT) then
          rcode = pio_get_att(pioid, varid, "scale_factor", scale_factor)
          if(rcode /= PIO_NOERR) then
             call shr_sys_abort('DATATYPE PIO_SHORT requires attributes scale_factor')
          endif
          rcode = pio_get_att(pioid, varid, "add_offset", add_offset)
          if(rcode /= PIO_NOERR) then
             call shr_sys_abort('DATATYPE PIO_SHORT requires attributes add_offset')
          endif
          rcode = pio_get_att(pioid, varid, "_FillValue", fillvalue_i2)
       endif
       if(rcode == PIO_NOERR) handlefill=.true.
       call PIO_seterrorhandling(pioid, old_error_handle)

       if (debug>0 .and. sdat%mainproc)  then
          write(sdat%stream(1)%logunit,F02)' reading '//&
               trim(per_stream%fldlist_stream(nf))//' into '//trim(per_stream%fldlist_model(nf)),&
               ' at time index: ',nt
       end if

       ! read the data
       call pio_setframe(pioid, varid, int(nt,kind=Pio_Offset_Kind))
       if (pio_iovartype == PIO_REAL) then
          ! -----------------------------
          ! pio_iovartype is PIO_REAL
          ! -----------------------------
          if (stream_nlev > 1) then
             if (per_stream%stream_pio_iodesc_set) then
                call pio_read_darray(pioid, varid, per_stream%stream_pio_iodesc, data_real2d, rcode)
             else
                rcode = pio_get_var(pioid, varid,start=(/1,1,1,nt/),count=(/1,1,1,1/), ival=data_real2d)
             end if
             if ( rcode /= PIO_NOERR ) then
                call shr_sys_abort(' ERROR: reading in variable: '// trim(per_stream%fldlist_stream(nf)))
             end if
             if (handlefill) then
                ! Single point streams are not allowed to have missing values
                if (stream%mapalgo == 'none' .and. any(data_real2d == fillvalue_r4)) then
                   write(errmsg,*) ' ERROR: _Fillvalue found in stream input variable: '// trim(per_stream%fldlist_stream(nf))
                   if(sdat%mainproc) write(sdat%stream(1)%logunit,*) trim(errmsg)
                   call shr_sys_abort(errmsg)
                endif
                do lev = 1,stream_nlev
                   do n = 1,size(dataptr2d, dim=2)
                      if (.not. shr_infnan_isnan(data_real2d(n,lev)) .and. data_real2d(n,lev) .ne. fillvalue_r4) then
                         dataptr2d(lev,n) = real(data_real2d(n,lev), kind=r8) ! Note the order of indices
                      else
                         dataptr2d(lev,n) = r8fill
                      endif
                   enddo
                end do
             else
                do lev = 1,stream_nlev
                   do n = 1,size(dataptr2d, dim=2)
                      dataptr2d(lev,n) = real(data_real2d(n,lev), kind=r8)
                   end do
                end do
             end if
          else ! stream_nlev == 1
             if (per_stream%stream_pio_iodesc_set) then
                call pio_read_darray(pioid, varid, per_stream%stream_pio_iodesc, data_real1d, rcode)
             else
                rcode = pio_get_var(pioid, varid,start=(/1,1,nt/),count=(/1,1,1/), ival=data_real1d)
             endif
             if ( rcode /= PIO_NOERR ) then
                call shr_sys_abort(' ERROR: reading in variable: '// trim(per_stream%fldlist_stream(nf)))
             end if
             if (handlefill) then
                ! Single point streams are not allowed to have missing values
                if (stream%mapalgo == 'none' .and. any(data_real1d == fillvalue_r4)) then
                   write(errmsg,*) ' ERROR: _Fillvalue found in stream input variable: '// trim(per_stream%fldlist_stream(nf))
                   if(sdat%mainproc) write(sdat%stream(1)%logunit,*) trim(errmsg)
                   call shr_sys_abort(errmsg)
                endif

                do n=1,size(dataptr1d)
                   if(.not. shr_infnan_isnan(data_real1d(n)) .and. data_real1d(n) .ne. fillvalue_r4) then
                      dataptr1d(n) = real(data_real1d(n), kind=r8)
                   else
                      dataptr1d(n) = r8fill
                   endif
                enddo
             else
                dataptr1d(:) = real(data_real1d(:),kind=r8)
             endif
          end if

       else if (pio_iovartype == PIO_DOUBLE) then
          ! -----------------------------
          ! pio_iovartype is PIO_DOUBLE
          ! -----------------------------
          if (stream_nlev > 1) then
             if (per_stream%stream_pio_iodesc_set) then
                call pio_read_darray(pioid, varid, per_stream%stream_pio_iodesc, data_dbl2d, rcode)
             else
                rcode = pio_get_var(pioid, varid,start=(/1,1,1,nt/), count=(/1,1,1,1/), ival=data_dbl2d)
             end if
             if ( rcode /= PIO_NOERR ) then
                call shr_sys_abort(' ERROR: reading in 2d double variable: '// trim(per_stream%fldlist_stream(nf)))
             end if
             if (handlefill) then
                ! Single point streams are not allowed to have missing values
                if (stream%mapalgo == 'none' .and. any(data_dbl2d == fillvalue_r8)) then
                   write(errmsg,*) ' ERROR: _Fillvalue found in stream input variable: '// trim(per_stream%fldlist_stream(nf))
                   if(sdat%mainproc) write(sdat%stream(1)%logunit,*) trim(errmsg)
                   call shr_sys_abort(errmsg)
                endif
                do lev = 1,stream_nlev
                   do n = 1,size(dataptr2d, dim=2)
                      if (.not. shr_infnan_isnan(data_dbl2d(n,lev)) .and. data_dbl2d(n,lev) .ne. fillvalue_r8) then
                         dataptr2d(lev,n) = data_dbl2d(n,lev)
                      else
                         dataptr2d(lev,n) = r8fill
                      endif
                   enddo
                end do
             else
                do lev = 1,stream_nlev
                   do n = 1,size(dataptr2d, dim=2)
                      dataptr2d(lev,n) = data_dbl2d(n,lev)
                   end do
                end do
             end if
          else ! stream_nlev == 1
             if (per_stream%stream_pio_iodesc_set) then
                call pio_read_darray(pioid, varid, per_stream%stream_pio_iodesc, data_dbl1d, rcode)
             else
                rcode = pio_get_var(pioid, varid,start=(/1,1,nt/), count=(/1,1,1/), ival=data_dbl1d)
             endif
             if ( rcode /= PIO_NOERR ) then
                call shr_sys_abort(' ERROR: reading in variable: '// trim(per_stream%fldlist_stream(nf)))
             end if
             if (handlefill) then
                ! Single point streams are not allowed to have missing values
                if (stream%mapalgo == 'none' .and. any(data_dbl1d == fillvalue_r8)) then
                   write(errmsg,*) ' ERROR: _Fillvalue found in stream input variable: '// trim(per_stream%fldlist_stream(nf))
                   if(sdat%mainproc) write(sdat%stream(1)%logunit,*) trim(errmsg)
                   call shr_sys_abort(errmsg)
                endif
                do n = 1,size(dataptr1d)
                   if (.not. shr_infnan_isnan(data_dbl1d(n)) .and. data_dbl1d(n) .ne. fillvalue_r8) then
                      dataptr1d(n) = data_dbl1d(n)
                   else
                      dataptr1d(n) = r8fill
                   end if
                enddo
             else
                do n = 1,size(dataptr1d)
                   dataptr1d(n) = data_dbl1d(n)
                end do
             endif
          end if

       elseif (pio_iovartype == PIO_SHORT) then
          ! -----------------------------
          ! pio_iovartype is PIO_SHORT
          ! -----------------------------
          if (stream_nlev > 1) then
             if (per_stream%stream_pio_iodesc_set) then
                call pio_read_darray(pioid, varid, per_stream%stream_pio_iodesc, data_short2d, rcode)
             else
                rcode = pio_get_var(pioid, varid,start=(/1,1,1,nt/), count=(/1,1,1,1/), ival=data_short2d)
             end if
             if ( rcode /= PIO_NOERR ) then
                call shr_sys_abort(' ERROR: reading in 2d short variable: '// trim(per_stream%fldlist_stream(nf)))
             end if
             if (handlefill) then
                do lev = 1,stream_nlev
                   do n = 1,lsize
                      if(data_short2d(n,lev) .eq. fillvalue_i2) then
                         dataptr2d(lev,n) = r8fill
                      else
                         dataptr2d(lev,n) = real(data_short2d(lev,n),r8) * scale_factor + add_offset
                      endif
                   enddo
                end do
             else
                do lev = 1,stream_nlev
                   do n = 1,lsize
                      dataptr2d(lev,n) = real(data_short2d(n,lev),r8) * scale_factor + add_offset
                   enddo
                end do
             end if
          else ! stream_nlev == 1
             if (per_stream%stream_pio_iodesc_set) then
                call pio_read_darray(pioid, varid, per_stream%stream_pio_iodesc, data_short1d, rcode)
             else
                rcode = pio_get_var(pioid, varid,start=(/1,1,nt/),count=(/1,1,1/), ival=data_short1d)
             endif
             if ( rcode /= PIO_NOERR ) then
                call shr_sys_abort(' ERROR: reading in variable: '// trim(per_stream%fldlist_stream(nf)))
             end if
             if (handlefill) then
                do n=1,lsize
                   if(data_short1d(n).eq.fillvalue_i2) then
                      dataptr1d(n) = r8fill
                   else
                      dataptr1d(n) = real(data_short1d(n),r8) * scale_factor + add_offset
                   endif
                enddo
             else
                do n=1,lsize
                   dataptr1d(n) = real(data_short1d(n),r8) * scale_factor + add_offset
                enddo
             endif
          end if

       else
          ! -----------------------------
          ! pio_iovartype is not supported
          ! -----------------------------
          call shr_sys_abort(subName//"ERROR: only double, real and short types are supported for stream read")

       end if

       if(associated(dataptr2d_src) .and. trim(per_stream%fldlist_model(nf)) .eq. uname) then
          ! save in dataptr2d_src
          dataptr2d_src(1,:) = dataptr1d(:)
       elseif(associated(dataptr2d_src) .and. trim(per_stream%fldlist_model(nf)) .eq. vname) then
          dataptr2d_src(2,:) = dataptr1d(:)
       else if (per_stream%stream_pio_iodesc_set) then
          ! Regrid the field_stream read in to the model mesh
          call dshr_fldbun_getfieldN(fldbun_data, nf, field_dst, rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
          call ESMF_FieldRegrid(per_stream%field_stream, field_dst, routehandle=per_stream%routehandle, &
               termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=.false., zeroregion=ESMF_REGION_TOTAL, rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       else
          call dshr_fldbun_getfieldN(fldbun_data, nf, field_dst, rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
          call ESMF_FieldFill(field_dst, dataFillScheme="const", const1=dataptr1d(1), rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       endif
    enddo

    ! Both components of a vector stream must be in the same input stream file
    if (associated(dataptr2d_src) .and. associated(dataptr1d)) then

       ! get lon and lat of stream u and v fields
       lsize = size(dataptr1d)
       allocate(dataptr(lsize))

       call ESMF_MeshGet(per_stream%stream_mesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       allocate(nu_coords(spatialDim*numOwnedElements))
       call ESMF_MeshGet(per_stream%stream_mesh, ownedElemCoords=nu_coords)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       allocate(nv_coords(spatialDim*numOwnedElements))
       call ESMF_MeshGet(per_stream%stream_mesh, ownedElemCoords=nv_coords)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       do i=1,lsize
          dataptr(i) = dataptr2d_src(1,i)
          lon = nu_coords(2*i-1)
          lat = nu_coords(2*i)
          sinlon = sin(lon*deg2rad); coslon = cos(lon*deg2rad)
          sinlat = sin(lat*deg2rad); coslat = cos(lat*deg2rad)
          dataptr2d_src(1,i) = (coslon * dataptr(i) - sinlon * dataptr2d_src(2,i))
          dataptr2d_src(2,i) = (sinlon * dataptr(i) + coslon * dataptr2d_src(2,i))
       enddo
       field_vector_dst = ESMF_FieldCreate(sdat%model_mesh, ESMF_TYPEKIND_r8, name='field_vector_dst', &
            ungriddedLbound=(/1/), ungriddedUbound=(/2/), gridToFieldMap=(/2/), meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       if (.not. ESMF_FieldIsCreated(per_stream%field_stream_vector)) then
          call shr_sys_abort('ERROR: per_stream%field_stream_vector has not been created')
       end if

       call ESMF_FieldRegrid(per_stream%field_stream_vector, field_vector_dst, per_stream%routehandle, &
            termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       call ESMF_FieldGet(field_vector_dst, farrayPtr=dataptr2d_dst, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       call dshr_fldbun_getFldPtr(fldbun_data, trim(uname), data_u_dst, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       call dshr_fldbun_getFldPtr(fldbun_data, trim(vname), data_v_dst, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       do i = 1,size(data_u_dst)
          lon = sdat%model_lon(i)
          lat = sdat%model_lat(i)
          sinlon = sin(lon*deg2rad); coslon = cos(lon*deg2rad)
          sinlat = sin(lat*deg2rad); coslat = cos(lat*deg2rad)
          data_u_dst(i) =  coslon * dataptr2d_dst(1,i) + sinlon * dataptr2d_dst(2,i)
          data_v_dst(i) = -sinlon * dataptr2d_dst(1,i) + coslon * dataptr2d_dst(2,i)
       enddo
       deallocate(dataptr)
    endif

    if (pio_iovartype == PIO_REAL) then
       if (allocated(data_real1d)) deallocate(data_real1d)
       if (allocated(data_real2d)) deallocate(data_real2d)
    else if (pio_iovartype == PIO_DOUBLE) then
       if (allocated(data_dbl1d)) deallocate(data_dbl1d)
       if (allocated(data_dbl2d)) deallocate(data_dbl2d)
    else if (pio_iovartype == PIO_SHORT) then
       if (allocated(data_short1d)) deallocate(data_short1d)
       if (allocated(data_short2d)) deallocate(data_short2d)
    endif
    if (.not. per_stream%stream_pio_iodesc_set) then
       deallocate(dataptr1d)
    endif
    call ESMF_TraceRegionExit(trim(istr)//'_readpio')

  end subroutine shr_strdata_readstrm

  !===============================================================================
  subroutine shr_strdata_set_stream_iodesc(sdat, per_stream, fldname, pioid, rc)

    ! Set stream_pio_iodesc and stream_pio_iodesc_set in stream

    ! input/output variables
    type(shr_strdata_type)      , intent(in)    :: sdat  ! strdata data data-type
    type(shr_strdata_perstream) , intent(inout) :: per_stream
    character(len=*)            , intent(in)    :: fldname
    type(file_desc_t)           , intent(inout) :: pioid
    integer                     , intent(out)   :: rc

    ! local variables
    integer                 :: stream_nlev
    integer                 :: gsize2d
    integer                 :: pio_iovartype
    integer                 :: n, m, cnt
    type(var_desc_t)        :: varid
    integer                 :: ndims
    character(CS)           :: dimname
    integer, allocatable    :: dimids(:)
    integer, allocatable    :: dimlens(:)
    type(ESMF_DistGrid)     :: distGrid
    integer                 :: lsize
    integer, pointer        :: compdof(:)
    integer, pointer        :: compdof3d(:)
    integer                 :: rCode ! pio return code (only used when pio error handling is PIO_BCAST_ERROR)
    character(*), parameter :: subname = '(shr_strdata_set_stream_iodesc) '
    character(*), parameter :: F00  = "('(shr_strdata_set_stream_iodesc) ',a,i8,2x,i8,2x,a)"
    character(*), parameter :: F01  = "('(shr_strdata_set_stream_iodesc) ',a,i8,2x,i8,2x,a)"
    character(*), parameter :: F02  = "('(shr_strdata_set_stream_iodesc) ',a,i8,2x,i8,2x,i8,2x,a)"
    !-------------------------------------------------------------------------------

    rc = ESMF_SUCCESS

    ! nullify local pointers
    nullify(compdof)
    nullify(compdof3d)

    ! set the number of vertical levels to a local variable
    stream_nlev = per_stream%stream_nlev

    ! query the variable fldname in the stream dataset for its dimensions
    rcode = pio_inq_varid(pioid, trim(fldname), varid)
    rcode = pio_inq_varndims(pioid, varid, ndims)

    ! allocate memory for dimids and dimlens
    allocate(dimids(ndims))
    allocate(dimlens(ndims))
    rcode = pio_inq_vardimid(pioid, varid, dimids(1:ndims))
    do n = 1, ndims
       rcode = pio_inq_dimlen(pioid, dimids(n), dimlens(n))
    end do

    ! determine compdof for stream
    call ESMF_MeshGet(per_stream%stream_mesh, elementdistGrid=distGrid, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    call ESMF_DistGridGet(distGrid, localDe=0, elementCount=lsize, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    allocate(compdof(lsize))
    call ESMF_DistGridGet(distGrid, localDe=0, seqIndexList=compdof, rc=rc)
    if (ChkErr(rc,__LINE__,u_FILE_u)) return
    if (stream_nlev > 1) then
       allocate(compdof3d(stream_nlev*lsize))
       ! Assume that first 2 dimensions correspond to the compdof
       gsize2d = dimlens(1)*dimlens(2)
       cnt = 0
       do n = 1,stream_nlev
          do m = 1,size(compdof)
             cnt = cnt + 1
             compdof3d(cnt) = (n-1)*gsize2d + compdof(m)
          enddo
       enddo
    end if

    ! determine type of the variable
    rcode = pio_inq_vartype(pioid, varid, pio_iovartype)

    ! determine io descriptor
    if (ndims == 2) then
       if (sdat%mainproc) then
          write(sdat%stream(1)%logunit,F00) 'setting iodesc for : '//trim(fldname)// &
               ' with dimlens(1), dimlens2 = ',dimlens(1),dimlens(2),&
               ' variable has no time dimension '
       end if
       call pio_initdecomp(sdat%pio_subsystem, pio_iovartype, (/dimlens(1),dimlens(2)/), compdof, &
            per_stream%stream_pio_iodesc)

    else if (ndims == 3) then
       rcode = pio_inq_dimname(pioid, dimids(ndims), dimname)
       if (stream_nlev > 1) then
          write(sdat%stream(1)%logunit,F01) 'setting iodesc for : '//trim(fldname)// &
               ' with dimlens(1), dimlens(2), dimlens(3) = ',dimlens(1),dimlens(2), dimlens(3), &
               ' variable has no time dimension '//trim(dimname)
          call pio_initdecomp(sdat%pio_subsystem, pio_iovartype, (/dimlens(1),dimlens(2),dimlens(3)/), compdof3d, &
               per_stream%stream_pio_iodesc)
       else if (trim(dimname) == 'time' .or. trim(dimname) == 'nt') then
          if (sdat%mainproc) then
             write(sdat%stream(1)%logunit,F01) 'setting iodesc for : '//trim(fldname)// &
                  ' with dimlens(1), dimlens(2) = ',dimlens(1),dimlens(2),&
                  ' variable as time dimension '//trim(dimname)
          end if
          call pio_initdecomp(sdat%pio_subsystem, pio_iovartype, (/dimlens(1),dimlens(2)/), compdof, &
               per_stream%stream_pio_iodesc)
       end if

    else if (ndims == 4) then
       rcode = pio_inq_dimname(pioid, dimids(ndims), dimname)
       if (stream_nlev > 1 .and. (trim(dimname) == 'time' .or. trim(dimname) == 'nt')) then
          if (sdat%mainproc) then
             write(sdat%stream(1)%logunit,F02) 'setting iodesc for : '//trim(fldname)// &
                  ' with dimlens(1), dimlens(2),dimlens(3) = ',dimlens(1),dimlens(2),dimlens(3),&
                  ' variable has time dimension '
          end if
          call pio_initdecomp(sdat%pio_subsystem, pio_iovartype, (/dimlens(1),dimlens(2),dimlens(3)/), compdof3d, &
               per_stream%stream_pio_iodesc)
       else
          write(6,*)'ERROR: dimlens= ',dimlens
          call shr_sys_abort(trim(subname)//' dimlens = 4 assumes a time dimension')
       end if

    else
       write(6,*)'ERROR: dimlens= ',dimlens
       call shr_sys_abort(trim(subname)//' only ndims of 2 and 3 and 4 are currently supported')
    end if

    ! deallocate memory
    deallocate(compdof)
    if (associated(compdof3d)) deallocate(compdof3d)
    deallocate(dimids)
    deallocate(dimlens)

  end subroutine shr_strdata_set_stream_iodesc

  !===============================================================================
  subroutine shr_strdata_get_stream_pointer_1d(sdat, strm_fld, strm_ptr, rc)

    ! Set a pointer, strm_ptr, for field, strm_fld, into sdat fldbun_model field bundle

    ! input/output variables
    type(shr_strdata_type) , intent(in)    :: sdat
    character(len=*)       , intent(in)    :: strm_fld
    real(r8)               , pointer       :: strm_ptr(:)
    integer                , intent(out)   :: rc

    ! local variables
    integer :: ns, nf
    logical :: found
    character(len=*), parameter :: subname='(shr_strdata_get_stream_pointer_1d)'
    character(*)    , parameter :: F00 = "('(shr_strdata_get_stream_pointer_1d) ',8a)"
    ! ----------------------------------------------

    rc = ESMF_SUCCESS

    ! loop over all input streams and determine if the strm_fld is in the field bundle of the target stream
    do ns = 1, shr_strdata_get_stream_count(sdat)
       found = .false.
       ! Check if requested stream field is read in - and if it is then point into the stream field bundle
       do nf = 1,size(sdat%pstrm(ns)%fldlist_model)
          if (trim(strm_fld) == trim(sdat%pstrm(ns)%fldlist_model(nf))) then
             call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_model, trim(sdat%pstrm(ns)%fldlist_model(nf)), &
                  fldptr1=strm_ptr, rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
             if (sdat%mainproc) then
                write(sdat%stream(1)%logunit,F00)' strm_ptr is allocated for stream field strm_'//trim(strm_fld)
             end if
             found = .true.
             exit
          end if
       end do
       if (found) exit
    end do
  end subroutine shr_strdata_get_stream_pointer_1d

  !===============================================================================
  subroutine shr_strdata_get_stream_pointer_2d(sdat, strm_fld, strm_ptr, rc)

    ! Set a pointer, strm_ptr, for field, strm_fld, into sdat fldbun_model field bundle

    ! input/output variables
    type(shr_strdata_type) , intent(in)    :: sdat
    character(len=*)       , intent(in)    :: strm_fld
    real(r8)               , pointer       :: strm_ptr(:,:)
    integer                , intent(out)   :: rc

    ! local variables
    integer :: ns, nf
    logical :: found
    character(len=*), parameter :: subname='(shr_strdata_get_stream_pointer_2d)'
    character(*)    , parameter :: F00 = "('(shr_strdata_get_stream_pointer_2d) ',8a)"
    ! ----------------------------------------------

    rc = ESMF_SUCCESS

    ! loop over all input streams and determine if the strm_fld is in the field bundle of the target stream
    do ns = 1, shr_strdata_get_stream_count(sdat)
       found = .false.
       ! Check if requested stream field is read in - and if it is then point into the stream field bundle
       do nf = 1,size(sdat%pstrm(ns)%fldlist_model)
          if (trim(strm_fld) == trim(sdat%pstrm(ns)%fldlist_model(nf))) then
             call dshr_fldbun_getfldptr(sdat%pstrm(ns)%fldbun_model, trim(sdat%pstrm(ns)%fldlist_model(nf)), &
                  fldptr2=strm_ptr, rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
             if (sdat%mainproc) then
                write(sdat%stream(1)%logunit,F00)' strm_ptr is allocated for stream field strm_'//trim(strm_fld)
             end if
             found = .true.
             exit
          end if
       end do
       if (found) exit
    end do
  end subroutine shr_strdata_get_stream_pointer_2d

end module dshr_strdata_mod
